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The  fact  that  'realistic'  galactic  potentials  are  not  necessarily  integrable 
implies  that  a  large  number  of  stars  should  be  expected  to  follow  stochastic 
trajectories.  Orbits  of  stars  in  an  elliptical  galaxy  are  assumed  to  be  affected 
by  both  the  smooth  mean  field  and  internal  fluctuations  associated  with  dis- 
creteness effects.  The  interplay  of  these  two  effects  is  studied  in  the  context  of 
the  gravitational  TV-body  problem.  These  studies  indicate  that  both  the  mean 
field  and  discreteness  effects  have  an  effect  on  the  instability  of  the  trajectories. 
These  discreteness  effects  are  also  modelled  as  noise,  perturbing  trajectories  in 
a  nonintegrable  potential.  The  conclusion  here  is  that  perturbations  in  the 
phase  space  variables  grow  exponentially  on  a  time  scale  associated  with  the 
instability  of  the  deterministic  trajectories,  even  though  the  time  scale  associ- 
ated with  changes  in  collisionless  invariants  remains  the  conventional  relaxation 
time.  Because  physically  relevant  nonintegrable  potentials  are  typically  non- 
hyperbolic,  the  standard  argument  that  shadowing  implies  that  noisy  orbits 
behave  like  deterministic  orbits  does  not  hold.  A  nonhyperbolic  shadowing 
theorem  is  used  to  determine  bounds  on  the  effect  that  noise  has  in  such  a 
potential.    This  noise  is  interpreted  to  be  due  to  discreteness  effects  caused 
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by  close  encounters.  Shadowing  times  are  calculated  for  both  maps  and  dif- 
ferential equations,  yielding  the  following  conclusions.  First,  shadowing  times 
for  such  systems  tend  towards  a  roughly  log-normal  distribution.  These  dis- 
tributions appear  to  exhibit  a  fractal  dependence  on  initial  conditions.  It  has 
been  conjectured  previously  that  the  shadowing  times  are  expected  to  be  pro- 
portional to  1/Vo,  where  8  is  the  amplitude  of  the  noise.  It  is  found  that, 
for  low  noise  levels  in  Hamiltonian  differential  equations,  the  shadowing  times 
are  considerably  shorter  than  this  conjecture  would  indicate.  Finally,  the  loca- 
tions where  shadowing  breaks  down  appear  to  be  correlated  with  the  location 
of  KAM  tori  and  cantori  of  Hamiltonian  systems. 


CHAPTER  1 
INTRODUCTION  TO  GALACTIC  DYNAMICS 

Galaxies  are  systems  which  consist  of  a  large  number  of  stars,  gas,  and 
dust.  Typically,  the  number  of  stars  ranges  from  1010  to  1012.  Despite  the 
complicated  interactions  among  the  various  components  of  the  innumerable 
galaxies,  they  can  generally  be  classified  into  spirals  and  ellipticals,  with  a  few 
being  irregular.  Elliptical  galaxies  are  observed  to  be  virtually  dust-free,  and 
hence  are  usually  modelled  as  a  collection  of  stars  interacting  via  gravitational 
forces,  described  by  Newton's  Law  of  Gravitation.  As  such,  it  is  tempting 
to  describe  the  dynamics  of  these  galactic  systems  using  results  from  celestial 
mechanics,  classical  statistical  mechanics,  and  plasma  physics.  However,  there 
are  difficulties  with  each  of  these  approaches,  and  hence  galactic  dynamics 
should  be  considered  a  branch  of  theoretical  physics  in  its  own  right. 

The  analogy  with  calculations  in  celestial  mechanics,  which  describe  the 
motions  of  objects  in  the  solar  system,  may  appear  useful  since  we  are  consid- 
ering gravitationally  bound  systems.  However,  they  are  based  on  perturbation 
expansions  which  are  based  on  a  hierarchy  of  mass  scales  that  do  not  converge 
in  the  case  of  stellar  systems.  Hence  any  mathematical  formalism  based  on 
celestial  mechanics  is  essentially  of  no  use  in  the  study  of  stellar  systems. 

The  large  number  of  particles  in  galactic  systems  requires  that  orbits  in 
a  galactic  system  be  treated  statistically.    However,  the  long-range  nature  of 
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gravitational  forces  means  that  they  must  be  treated  differently  than  the  short- 
range  forces  associated,  e.g.,  with  the  motions  of  gas  molecules  confined  to  a 
box.  In  a  later  chapter,  we  will  see  in  more  detail  why  orbits  in  a  gravitational 
system  do  not  obey  Boltzmann's  ergodic  hypothesis,  which  is  the  foundation 
of  equilibrium  statistical  mechanics. 

Finally,  plasma  physics,  which  also  involves  the  study  of  a  large  number 
of  particles  interacting  via  long-range  forces,  in  this  case  electrostatic,  requires 
many  of  the  same  mathematical  tools  as  galactic  dynamics.  However,  the  ex- 
istence of  both  positive  and  negative  charges  leads  to  screened  potentials  and 
static,  homogeneous  equilibria,  which  are  absent  in  galactic  systems.  Never- 
theless, the  collisional  Boltzmann  equation,  primarily  used  in  plasma  physics, 
is  frequently  relied  upon  in  the  study  of  galactic  systems  and,  in  fact,  Jeans 
formulated  the  equation  in  the  context  of  gravitational  systems  before  it  was 
used  for  plasma  physics. 

The  dynamics  of  stars  in  galactic  systems  may  be  extremely  complicated, 
and  any  reasonable  model  of  a  galactic  potential  must  take  into  account  the  fact 
that  individual  stellar  trajectories  may  not,  in  generic  situations,  be  calculated 
analytically.  Nevertheless,  methods  exist  that  make  a  statistical  treatment  of 
stellar  dynamics  feasible.  Standard  methods  include  N-body  calculations  and 
the  study  of  the  collisionless  Boltzmann  equation.  While  both  of  these  methods 
have  their  shortcomings,  to  be  described  later,  they  do  have  the  ability  to  yield 
useful  statistical  results. 

It  is  the  purpose  of  this  chapter  to  outline  various  results  of  galactic  dynam- 
ics, found  via  standard  methods,  and  to  motivate  the  idea  that  the  modern 
techniques  of  dynamical  systems  theory  are  necessary.    To  begin  with,  the 
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standard  paradigm  of  decoupling  collisional  effects  from  bulk  mean  field  effects 
will  be  discussed,  leading  to  the  conventional  calculation  of  a  characteristic 
relaxation  time  for  stellar  systems.  Then,  some  discussion  of  how  collisional 
effects  should  be  incorporated  will  be  included.  After  this,  I  will  attempt  to 
convince  the  reader  that  the  apparent  existence  of  nonintegrable  galactic  poten- 
tials makes  modern  dynamical  systems  techniques  of  use  for  galactic  systems. 
It  should  be  noted  that  the  arguments  presented  in  this  chapter  are  not 
rigorous.  The  field  of  galactic  dynamics  draws  upon  results  from  other  dis- 
ciplines, for  example,  plasma  physics,  which  are  based  on  assumptions  which 
may  not  be  suitable  for  galactic  systems.  Therefore,  the  bulk  of  this  chapter 
should  be  taken  as  an  attempt  to  briefly  discuss  the  conventional  wisdom  in 
the  field,  and  is  not  meant  to  imply  that  the  author  concurs  with  the  beliefs 
presented  here. 

1.1  Characteristic  Time  Scales  in  Stellar  Dynamics 

1.1.1  The  Gravothermal  Catastrophe 

When  one  desires  to  calculate  the  dynamical  properties  of  a  given  physical 
system  statistically,  it  is  usually  necessary  to  determine  the  proper  time  scales 
over  which  the  system  evolves.  One  is  often  interested,  for  example,  in  the 
approach  to  equilibrium.  However,  for  self-gravitating  systems,  care  must  be 
taken  in  defining  'equilibrium,'  since  a  thermal  equilibrium  does  not  exist  for 
such  systems. 

This  is  easily  seen,  as  discussed  in  detail,  e.g.,  by  Lynden-Bell  and  Lynden- 
Bell  [1].  Let  us  consider  the  following  heuristic  argument.  If  the  stars  in 
a  galaxy  were  considered  to  be  an  ideal  gas,  the  theorem  of  equipartition  of 
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energy  would  be  expected  to  hold,  relating  the  mean-square  velocity  of  the 
particles  v2  to  the  temperature  of  the  system  T: 

\m^=\kBT,  (1.1) 

where  m  is  the  stellar  mass  and  ks  is  Boltzmann's  constant.  The  Virial  Theo- 
rem, which  states  that  for  an  isolated  system,  neglecting  boundaries,  the  total 
energy  is  equal  to  the  negative  of  the  total  kinetic  energy,  allows  us  to  write  the 
total  energy  of  a  galactic  system  with  N  particles  in  terms  of  the  temperature 
of  the  system, 

E  =  ~NkBT.  (1.2) 

The  heat  capacity  of  the  system  is  then  calculated  to  be  negative: 

Cv^%  =  -\NkB-  (L3) 

This  apparently  paradoxical  result  of  a  negative  heat  capacity  is  a  conse- 
quence of  the  fact  that,  in  general,  gravitational  systems  cannot  be  described 
in  terms  of  a  canonical  ensemble,  and  cannot  be  considered  to  be  in  a  ther- 
mal equilibrium.  This  is  true  for  any  bound,  finite,  unconstrained,  gravitating 
system,  and  is  responsible  for  the  so-called  gravothermal  catastrophe:  a  ther- 
mally isolated  sphere  is  unstable  to  perturbations.  If  the  central  core  of  the 
sphere  becomes  hotter  than  the  periphery,  it  will  continue  to  become  hotter 
and  denser,  leading  to  a  runaway  effect.  Thus  one  cannot,  in  a  true  sense,  speak 
of  a  thermal  equilibrium  for  gravitational  systems.  Instead,  one  must  be  con- 
tent to  consider  gravitational  systems  that  approach,  for  example,  a  dynamical 
equilibrium.  It  is  this  sort  of  'quasi-equilibrium'  that  will  be  considered  when 
we  wish  to  discuss  the  relevant  time  scales  of  the  gravitational  system. 
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1.1.2  The  Dynamical  Time 

Now  that  we  see  that  the  concept  of  'equilibrium'  is  an  elusive  quantity,  we 
must  be  careful  with  determining  the  time  scales  associated  with  an  approach  to 
equilibrium.  Let  us  consider  that  the  system  approaches  some  'quasi-stationary 
state'  due  to  the  action  of  the  long-range  gravitational  forces  acting  as  a  smooth 
mean  field.  The  system  can  then  be  considered  to  be  in  a  virial  equilibrium, 
so  the  kinetic  energy  of  the  system  should  be  comparable  in  magnitude  to  the 
potential  energy,  and 


v    ~ 


GNm 


(1.4) 


R 

Since  the  'collective'  force  acting  on  a  test  particle  of  unit  mass  will  have  a  typ- 
ical magnitude  of  Ff0t  ~  G^m,  we  can  determine  this  characteristic  dynamic 

time  scale  to  be  tj  ~  =L-,  which  is  the  crossing  time,  rc  ~  R/v,  where  R  is  the 

Ftot 

characteristic  length  scale  of  the  system  and  iJ  is  the  stellar  root-mean-square 
velocity. 

1.1.3  The  Binary  Relaxation  Time 

Of  course,  one  has  to  realize  that  not  only  the  bulk  mean  field,  but  also 
encounters  with  nearby  stars  can  affect  the  stellar  trajectories.  Generally,  one 
may  wish  to  calculate  the  'relaxation  time,'  the  characteristic  time  for  which 
a  star  would  forget  its  initial  conditions  via  interactions  with  nearby  stars. 
Specifically,  one  can  calculate  the  binary  relaxation  time,  a.k.a.  gravitational 
Rutherford  scattering,  which  is  usually  defined  to  be  the  time  required  for 
the  velocity  of  a  test  star  to  change  appreciably  due  to  graininess  effects  as  it 
crosses  the  system,  or 

rr  =  £j-rc,  (1.5) 
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where  Av±  is  the  integrated  deflection  in  the  velocity  due  to  the  rest  of  the 
N  stars.  If  we  again  utilize  the  virial  theorem,  it  can  be  shown  that  (see,  e.g. 

[2]) 

rr  wO.liV/lnA,  (1.6) 

where  A  is  the  ratio  of  the  impact  parameters  f"^*-.  In  plasma  physics,  bmax  is 
the  Debye  cutoff,  which  does  not  exist  in  gravitational  physics.  Instead,  bmax 
is  considered  to  be  about  the  size  of  the  system.  Exactly  how  one  defines  this 
value  is  open  to  debate,  however,  in  any  event  In  A  ~  In  TV.  The  determina- 
tion of  this  binary  relaxation  time,  however,  is  not  unique.  rr  may  also  be 
interpreted  as  the  time  over  which  appreciable  changes  in  energy  occur,  for  ex- 
ample, or  the  time  required  for  equipartition  of  energy  to  occur.  The  standard 
analysis  neglects  the  movement  of  the  other  stars  as  well.  In  any  event,  if  the 
masses  of  the  stars  are  all  comparable,  these  differences  do  not  have  much  of 
an  effect  on  the  value  of  rr .  A  typical  value  of  rr  for  an  elliptical  galaxy  would 
be  1014-1016  years. 

Numerical  values  for  these  characteristic  times  are  easy  to  determine.  For 
a  typical  elliptical  galaxy  of  N  W  1011  stars,  the  dynamical  time  is  of  the  order 
108  years,  whereas  the  binary  relaxation  time,  1014  years,  is  much  longer  than 
the  age  of  the  universe.  Thus,  from  this  analysis,  it  would  appear  that,  over 
relevant  time  scales,  binary  encounters  are  unimportant.  However,  it  must 
be  noted  that  this  derivation  of  the  binary  relaxation  time  is  based  on  the 
assumption  that,  without  encounters,  the  stellar  trajectories  follow  a  smooth, 
rectilinear  trajectory,  i.e.  the  trajectories  are  not  exponentially  unstable. 
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1.1.4  Collective  Relaxation 

Perhaps  a  more  relevant  concept  is  the  idea  of  collective  relaxation.  Here, 
we  are  assuming  that  the  system  exhibits  sensitive  dependence  on  initial  con- 
ditions, which  is  certainly  true  for  a  large  choice  of  potentials.  In  general,  this 
sensitivity  is  exponential,  causing  trajectories  initially  a  distance  x(0)  apart 
to  be,  after  a  time  i,  a  distance  e'/r*x(0)  apart.  This  time  scale,  the  col- 
lective relaxation  time  scale,  can  be  estimated  numerically  (see  Chapter  3) 
or  analytically.  The  mathematical  formulation  of  the  problem  was  discussed 
in  a  gravitational  setting  by  Kandrup  [3],  who  viewed  the  dynamics  of  the 
stellar  trajectories  as  a  geodesic  flow  on  a  3iV-dimensional  manifold  with  aver- 
age negative  curvature,  and  concluded  that  this  collective  relaxation  should  be 
manifest  on  a  time  scale  ~  rc.  This  conclusion  was  also  reached  by  Goodman, 
Heggie,  and  Hut  [4],  using  a  very  different  approach. 

It  should  be  stressed  again  at  this  point  that  dynamically  galaxies  are  rel- 
atively young  objects,  having  lived  for  typically  no  more  than  100  crossing 
times,  and  not  nearly  old  enough  for  binary  relaxation  mechanisms  to  have 
had  an  effect.  This  latter  point  is  especially  meaningful  when  one  considers 
the  effect  that  nearby  stars  may  have  on  the  divergence  of  trajectories  in  a  stel- 
lar system.  It  may  be  true  that  the  classical  binary  relaxation  time  proceeds 
on  a  time  scale  much  too  long  to  be  considered  for  galactic  systems.  However, 
this  does  not  a  priori  imply  that  nearby  encounters  are  irrelevant,  if  one  ap- 
propriately considers  the  collective  relaxation  of  the  system.  Nevertheless,  it 
is  commonly  assumed  that,  on  appropriate  time  scales,  nearby  stellar  encoun- 
ters are  unimportant,  an  assumption  that  would  make  analytical  calculations 
easier.  This  is  the  idea  explored  in  the  next  section. 


8 
1.2  The  Collisionless  Boltzmann  Equation 
It  is  common  to  assume  that,  if  a  gravitational  system  is  large,  the  effect 
of  close  encounters  between  nearby  stars  will  be  small  compared  to  the  overall 
mean  field.  That  this  be  so  can  presumably  be  determined  from  the  following 
argument.  Let  us  assume  here  that  all  stars  have  the  same  mass.  The  average 
distance  between  stars  in  a  system  of  size  R  is  R/N1'3,  yielding  a  potential 
energy  between  the  two  stars  E2  ^  GmNllz/R  per  star.  Thus,  the  ratio 
between  the  average  potential  energy  between  two  stars  and  the  total  potential 
energy  per  star,  Ej~  =  GmN/R,  is 

e2/et  a  i\r2/3. 

From  this,  it  is  frequently  believed  that  close  encounters  have  a  negligible  effect 
on  the  relaxation  of  galactic  systems. 

With  this  in  mind,  and  realizing  that  gravitational  systems  are  far  from 
being  in  thermodynamic  equilibrium,  it  is  standard  to  describe  the  distribution 
of  stars  in  phase  space  using  kinetic  equations.  The  state  of  the  system  can 
be  completely  described  by  the  iV-particle  distribution  function  /'  ',  which 
obeys  Liouville's  equation, 

¥-* 

where  -^  is  the  convective  derivative  in  the  iV-dimensional  phase  space.  The 
probability  that  a  star  will  be  found  in  a  certain  cell  in  phase  space  is  given  by 

f{N)(*l,. . . , xN,  pi, . . . ,  pN,  t)d3Xl  . . .  dZxNd3Pl  . . .  d3pN.  (1.8) 

It  should  be  apparent  that  pN)  is  of  little  use  when  N  is  large,  since  it  is 
such  a  complex  and  unwieldy  object.    Instead  of  working  directly  with  f^N\ 
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we  can  define  a  one-particle  distribution  function  by 

/W(*l,pi,i)  =  J  f{N\xi,. . . , xN,  Pl, . . . , pN, t)d3x2  . .  .  d3XNd3p2  . . .  d3pN, 

(1.9) 
where 

fW(x,p,t)d3xd3p  (1.10) 

is  the  probability  of  finding  a  single  particle  in  the  phase  space  volume  d3xd3p 
around  the  point  (x,  p). 

We  can  as  well  define  a  two-particle  distribution  function  as 

/(2)(xi,pi,X2,P2,0  =  /(1)(xi,Pl,0/(1)(x2,P2,<)[1  +  </(xl5Pl,X2,P2,0]5 

(1.11) 

where  g  is  a  two-particle  correlation  function.  If  we  then  make  the  claim  that 
g  =  0,  which  is  exact  for  noninteracting  particles  and  a  good  approximation  for 
weakly  interacting  particles,  we  obtain  for  the  one-particle  distribution  function 
(called  /  for  brevity) 

f  =  f +P.V/-V*.^  =  0.  (1.12) 

dt       at      m  dp 

Equation  (1.12)  is  called  the  collisionless  Boltzmann  equation,  which  is 
also  known  as  the  electrostatic  Vlasov  equation  in  plasma  physics.  In  stellar 
dynamics,  we  are  considering  only  the  gravitational  force  in  —  V$  and  the 
momentum  p,  whereas,  in  general,  electric  and  magnetic  forces  are  included  as 
well.  It  should  also  be  noted  that  the  gravitational  potential  $  depends  on  the 
mean  field,  and  can  be  self-consistently  determined  by  the  Poisson  equation 


V2$  =  -47rGm 


y"rfv/(x,v,t).  (1.13) 
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Unlike  the  full  N-dimensional  Liouville  equation,  the  collisionless  Boltz- 
raann  equation  is  simple  enough  to  be  tractable.  One  equilibrium  solution  of 
this  equation  is  simply  a  Maxwell- Boltzmann  distribution,  never  fully  realized 
for  gravitational  systems  due,  to  some  extent,  to  escapes  and  to  the  gravother- 
mal  catastrophe,  which  nevertheless  eventually  relax  to  a  quasi-Maxwellian 
equilibrium.  As  it  has  been  supposedly  determined  that  this  relaxation  could 
not  be  due  to  close  encounters,  another  relaxation  process  must  have  taken 
place  for  this  to  occur,  since  galaxies  are  apparently  observed  to  have  relaxed. 
This  collective  relaxation  has  been  called  'violent  relaxation'  by  Lynden-Bell 
[5],  for  which  the  distribution  function  /  changes  rapidly,  on  the  order  of  a 
dynamical  time  tj. 

The  existence  of  a  violent  relaxation,  however,  does  not  guarantee  a  ther- 
modynamic quasiequilibrium.  During  the  course  of  the  relaxation,  the  distri- 
bution of  particles  gets  scrambled.  It  is  then  useful  to  consider  a  coarse-grained 
distribution  function,  which  averages  over  phase  cells.  This  coarse-grained  dis- 
tribution function  then  obeys  an  exclusion  principle  akin  to  Fermi-Dirac  statis- 
tics, since  Liouville's  theorem  disallows  degeneracy  of  phase  elements.  How- 
ever, this  new  distribution  function  does  not  obey  the  collisionless  Boltzmann 
equation,  and  hence  its  equilibrium  solution  is  not  necessarily  Maxwellian. 

1.3  Collisional  Stellar  Dynamics 

1.3.1  The  Collisional  Boltzmann  Equation 

Of  course  the  collisionless  Boltzmann  equation  is  not  valid  whenever  close 
encounters  are  important.  In  this  case,  it  is  necessary  to  include  a  'collision 
term'  in  equation  (1.12),  noting  that  'collision'  does  not  imply  any  actual 
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physical  contact,  rather,  that  the  stars  are  close  enough  for  graininess  to  be 
relevant.  Thus  the  (collisional)  Boltzmann  equation  reads 

The  right-hand  side  of  equation  (1.14)  includes  only  instantaneous  binary 
collisions,  which  is  a  good  approximation  for  dilute  gasses  with  short  range 
forces,  and  plasmas,  for  which  the  long  range  interaction  is  screened.  Though 
the  assumptions  are  suspect,  it  is  still  generally  believed  to  be  relevant  for 
gravitational  systems  when  N  is  large. 

The  collision  term  is  nonzero  whenever  the  two-particle  direct  correlation 
function  is  nonzero;  loosely  speaking,  this  is  the  more-or-less  obvious  statement 
that  correlations  between  two  particles  are  due  to  a  close  encounter  between 
them.  When  this  is  true,  the  collision  term  in  the  Boltzmann  equation,  which 
describes  the  evolution  of  the  one-particle  distribution  function,  is  written  in 
terms  of  the  two-particle  distribution  function,  and,  in  general,  the  equation  for 
the  evolution  of  the  n-particle  distribution  function  contains  an  n  +  1-particle 
distribution  function.  This  hierarchy  of  equations  is  known  as  the  BBGKY 
hierarchy.  The  Boltzmann  equation,  the  first  of  the  equations  in  the  hierarchy, 
is  obtained  in  the  limit  of  dilute  gasses,  when  only  binary  collisions  need  be 
accounted  for. 

The  first  two  equations  of  the  BBGKY  hierarchy  are  relevant  and  useful 
for  plasma  physics,  due  to  the  short-range  nature  of  the  correlations  from  the 
Debye  screening.  However,  there  is  a  difficulty  with  using  these  equations  for 
gravitational  systems.  Due  to  the  long-range  nature  of  the  force,  the  mean 
field  can  never  be  neglected. 
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The  question  remains  how  to  treat  the  collision  term  in  the  Boltzmann 
equation.  This  term  can  be  approximated  if  we  consider  only  weak  encounters, 
that  is,  encounters  for  which  the  velocity  change  due  to  the  encounter  is  small 
compared  to  the  original  velocity  of  the  star.  This  is  a  reasonable  approxi- 
mation for  galactic  systems,  since  for  such  systems  encounters  which  yield  a 
drastic  velocity  change  are  extremely  rare.  We  can  then  expand  the  collision 
term  in  a  Taylor  series  for  the  phase  space  variables  of  the  star.  Truncating 
the  series  after  the  second-order  terms,  we  have  an  approximate  solution  for 
the  collision  term 

(f  l, = - 1  >« + \  t  iw'^M, 

1=1  i,j=i       j 

(1.15) 
where  w{  are  the  phase  space  variables. 

Equation  (1-15)  is  the  Fokker-Planck  approximation,  and  D(Awi)  and 
D(AwiAwj),  which  are  functional  of  /,  are  known  as  the  diffusion  coeffi- 
cients. The  difficulty  now  lies  in  calculating  the  diffusion  coefficients  from 
the  one-particle  distribution  function.  The  equation  is  similar  to  the  actual 
Fokker-Planck  equation,  a  master  equation  for  the  probability  P(y,  t) 

1.3.2  The  Langevin  Equation 

Collisional  effects  can  be  treated  with  an  entirely  different  approach  than 
through  the  Boltzmann  equation  as  considered  heretofore.  Consider  the  veloc- 
ity of  a  free  particle  (i.e.  in  the  absence  of  external  forces)  subject  to  random 
fluctuations.  The  equation  of  motion  of  the  particle  can  then  be  written 

v  =  -7v +  £(*).  (1.17) 
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This  is  the  classical  picture  of  Brownian  motion,  and  equation  (1.17)  is  known 
as  the  Langevin  equation.  The  force  on  the  right-hand  side  of  the  equation 
consists  of  a  drag  force  which  is  linear  in  the  velocity,  and  a  random  external 
force,  which  is  due  to  the  collisions  of  the  particle  with  other  particles.  L(t)  is 
considered  to  be  a  stochastic  process,  defined  so  that  the  average  (L(t))  =  0. 
Usually  when  dealing  with  rapidly  varying  stochastic  forces,  i.e.  the  random 
bombardment  of  a  Brownian  particle  by  the  atoms  in  the  surrounding  fluid, 
the  autocorrelation  function  is  expressed  as  a  delta  function, 

(L(t)L(t'))  =  T6(t  -  «'),  (1.18) 

where  T  is  a  constant. 

By  writing  equation  (1.17),  we  have  made  the  assumption  that  the  force 
acting  on  a  particle  can  be  decomposed  into  two  independent  parts,  one  be- 
ing a  smooth  friction,  and  the  other  due  to  the  discontinuities  of  the  matter 
distribution.  Two  assumptions  have  been  made  about  the  fluctuations:  L(t) 
is  assumed  to  be  independent  of  the  velocity  v,  and  L(t)  is  assumed  to  vary 
extremely  rapidly  compared  to  variations  in  v. 

The  above  assumptions  are  certainly  valid  for  the  case  of  Brownian  motion. 
In  stellar  dynamics,  one  may  estimate  the  the  average  period  of  the  fluctuation 
of  the  influence  of  the  immediate  neighborhood  of  stars  to  be  T  ~  D/v,  where 
t;  is  the  root  mean  square  relative  velocity  between  the  two  stars,  and  D  is 
the  average  distance  between  nearest  neighbors,  which  can  be  calculated  to  be 
D  ~  0.55396n-1'3,  where  n  is  the  average  number  of  stars  per  unit  volume. 
If  we  consider  the  effect  of  such  fluctuations  on  the  sun,  for  example,  we  have 
D  ~  3  parsecs,  and  v  ~  50  km/sec.    Thus  we  obtain  T  ~  6  x  104  years. 
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Comparing  this  to  the  period  of  galactic  rotation,  which  is  about  2  x  108  years, 
we  see  that  the  fluctuating  force  does  indeed  vary  sufficiently  rapidly. 

The  postulates  of  the  Langevin  equation  only  specify  the  first  two  moments, 
saying  nothing  of  the  higher  moments  of  L(t).  If  we  include  an  additional 
postulate  that  L(t)  is  Gaussian,  i.e.  all  higher  cumulants  are  zero,  L(t)  is 
called  Gaussian  white  noise.  The  solution  of  v(t)  can  then  be  shown  to  be 
Gaussian.  One  can  then  conclude  that  the  Langevin  equation  with  Gaussian 
white  noise,  determined  be  equation  (1.17)  and  equation  (1.18),  is  equivalent 
to  the  Fokker-Planck  equation 

dP(v,t)       d   _    rd2p 

where  7  and  F  are  related  via  a  fluctuation-dissipation  theorem. 

Equation  (1.17)  neglects  the  effects  of  the  mean  field,  however.  If  a  mean 
field  potential  V(r,  t)  is  included,  the  generalized  Langevin  equation  can  be 
written 

v  =  -7v  +  V(p,  t)  +  L(t).  (1.20) 

The  fact  that  this  mean  field  may  be  unstable  has  serious  implications  for  the 
relaxation  time  of  the  system,  which  will  be  discussed  in  Chapter  4. 

1.4  Evidence  for  Stochastic  Orbits  in  Galactic  Systems 
As  mentioned  previously,  elliptical  galaxies  have  classical  binary  relaxation 
time-scales  much  longer  than  the  age  of  the  universe.  However,  many  such 
galaxies  are  well  fit  by  isothermal  density  distributions,  as  if  they  were  fully 
relaxed.  Thus  it  is  apparent  that  the  classical  estimate  of  the  relaxation  time 
is  inadequate. 
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A  primary  assumption  of  essentially  all  classical  estimates  of  the  relax- 
ation time  is  the  existence  of  an  integrable  mean  field;  stars  follow  a  rectilinear 
trajectory  except  when  they  experience  a  two-body  encounter.  In  three  di- 
mensions, most  potentials  that  can  be  constructed  are  nonintegrable,  leading 
to  chaotic  behavior,  with  integrable  potentials  being  a  set  of  measure  zero. 
Unlike  integrable  systems,  which  contain  only  regular  orbits,  nonintegrable 
systems  generally  exhibit  both  regular  and  stochastic  orbits,  which  would  be 
expected  to  have  drastic  effects  on  the  dynamical  evolution  of  the  system.  A 
more  detailed  description  of  nonintegrable  systems  is  forthcoming  in  chapter 
2. 

Nonetheless,  it  is  still  standard  practice  to  study  potentials  which  admit 
only  regular  orbits,  the  class  of  Stackel  potentials  being  an  example  of  com- 
monly studied  models  of  triaxial  galaxies.  It  could  be  argued  that,  even  if  the 
potential  is  not  completely  integrable,  the  fraction  of  irregular  orbits  may  yet 
be  small.  However,  recent  evidence  suggests  that  at  least  a  reasonable  frac- 
tion of  galaxies  have  a  significant  fraction  of  stochastic  orbits.  Such  evidence 
includes  numerical  studies  of  'realistic'  galactic  potentials,  as  well  as  recent 
observational  evidence  obtained  from  the  Hubble  Space  Telescope. 

Spiral  galaxies  that  are  nearly  axisymmetric  may  be  argued  to  contain 
predominantly  regular  orbits.  However,  spiral  galaxies  are  frequently  observed 
that  contain  a  large  amplitude  central  bar.  It  has  been  shown  numerically  [6,7] 
that  the  existence  of  such  bars  dramatically  increases  the  fraction  of  stochastic 
orbits.  Models  of  barred  galaxies  which  also  include  a  central  mass  concen- 
tration were  studied  by  Hasan  and  Norman     [8].    They  conclude  that,  even 
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though  the  central  mass  tends  to  dissolve  the  bar,  orbits  close  to  the  center  of 
the  galaxy  tend  to  become  stochastic. 

Significant  deviations  from  integrability  are  found  in  models  of  elliptical 
galaxies  as  well.  Considering  the  simplicity  of  the  shape  of  elliptical  galaxies, 
one  may  be  inclined  to  expect  that  regular  orbits  are  prevalent.  This  is  not 
necessarily  the  case,  however,  and  while  it  often  would  seem  convenient  to 
construct  stationary  configurations  using  nothing  but  regular  orbits,  the  pres- 
ence of  stochastic  orbits  is  not  inconsistent  with  the  existence  of  stationary 
quasi-equilibria. 

As  has  long  been  known,  typical  galaxy  components  tend  to  be  triaxial 
[9,10],  and  perturbations  of  even  integrable  models  typically  lead  to  a  finite 
fraction  of  stochastic  orbits,  as  is  the  case  with  Stackel  potentials  perturbed 
by  a  small  central  singularity  [11].  A  substantial  increase  in  stochasticity  was 
also  observed  in  numerical  experiments  by  Udry  and  Pfenniger  [12]  with  the 
addition  of  a  perturbation.  In  these  studies,  triaxial  galaxies  included  either  a 
central  mass  perturbation,  or  a  multipolar  perturbation  which  may  model  the 
twisting  of  isophotes  that  is  frequently  observed. 

These  numerical  observations  are  backed  up  by  recent  results  from  the 
Hubble  Space  Telescope  (HST).  Observational  evidence  obtained  from  the  HST 
confirms  and  extends  the  results  of  ground-based  observations  of  giant  ellipti- 
cals, which  are  shown  to  be  nonrotating,  anisotropic,  and  moderately  triaxial, 
and  typically  to  have  cuspy  cores  [13].  These  cusps  are  due  to  large  central 
mass  concentrations,  possibly  black  holes.  Galaxies  which  are  observed  by  the 
HST  to  have  cusps  have  been  classified  as  either  'core'  galaxies  or  'power-law' 
galaxies  [14].  Galactic  models  with  similar  cusps  have  been  studied  by  Merritt 
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and  Fridman  [15].  These  authors  studied  strongly  triaxial  systems  with  weak 
cusps,  for  which  the  central  density  varied  as  p  ~  r  ,  and  strong  cusps,  with 
p  ~  r-2,  which  have  density  profiles  similar  to  those  of  the  'core'  and  'power- 
law'  galaxies  observed  by  HST.  The  results  indicate  that  about  80  percent  of 
the  orbits  with  zero  initial  velocity  are  stochastic  for  a  wide  variety  of  energies 
for  the  strong  cusp  model,  and,  surprisingly,  between  60  and  80  percent  of  the 
orbits  are  stochastic  even  for  the  weak  cusp  model,  depending  on  the  orbital 
energy.  While  these  results  are  for  systems  which  are  more  strongly  triaxial 
than  observed,  they  do  indicate  that  a  substantial  amount  of  orbits  in  a  typical 
realistic  giant  elliptical  galaxy  would  be  expected  to  be  stochastic. 

The  above  evidence  seems  to  indicate  that,  if  one  wishes  to  consider  orbits 
of  stars  in  a  realistic  galactic  potential,  dynamical  effects  due  to  a  finite  degree 
of  stochasticity  should  not  be  ignored.  In  the  next  chapter,  we  will  introduce 
concepts  from  dynamical  systems  theory  that  will  allow  us  to  adequately  study 
the  dynamics  of  particles  in  such  nonintegrable  systems. 


CHAPTER  2 

CHAOS  IN  HAMILTONIAN 

DYNAMICAL  SYSTEMS 

The  purpose  of  this  chapter  is  to  introduce  the  concepts  and  techniques  of 
dynamical  systems  theory  necessary  for  the  understanding  of  chaos  in  Hamil- 
tonian  systems.  Dynamical  systems  can  be  classified  as  either  continuous  sys- 
tems, usually  described  as  a  system  of  differential  equations,  or  discrete  systems 
(maps).  A  system  of  N  first-order,  autonomous,  ordinary  differential  equations 
can  be  written  as  dx/dt  =  .F[x(£)],  where  x  is  an  N-dimensional  vector.  Maps 
can  be  written  in  vector  form  as  xn+i  =  M(xn),  where  xn  has  N  components. 
Maps  are  generally  easier  to  study  numerically;  however,  there  is  no  universal 
way  to  convert  continuous  systems  into  discrete  systems  while  being  certain 
not  to  alter  the  dynamics  in  a  qualitatively  significant  way.  Nevertheless,  under 
certain  circumstances,  Poincare  surfaces  provide  a  useful  method  of  reducing 
an  N-dimensional  continuous  system  to  an  (N-l)- dimensional  discrete  system. 
These  will  be  treated  later  in  this  chapter. 

Let  us  restrict  our  attention  to  conservative,  i.e.  Hamiltonian,  systems, 
for  which  the  phase  volume  Vx  •  F(x)  =  0,  as  required  by  Liouville's  theorem. 
The  N-dimensional  space  through  which  the  system  evolves  is  called  phase 
space.  Given  an  initial  condition  x(0)  (or  xo  for  maps),  the  path  through  phase 
space  along  which  a  particle  will  travel  is  called  an  orbit  or  trajectory.  These 
terms  have  slightly  different  meanings  in  some  of  the  mathematical  literature; 
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however,  this  distinction  is  not  universal  and  from  here  on  the  terms  will  be 
used  interchangeably. 

2.1  Integrable  Hamiltonian  Dynamical  Systems 
A  Hamiltonian  mechanical  system  is  given  by  the  phase  space,  a  symplec- 
tic  structure  on  the  phase  space  given  by  the  symplectic  two- form  u2,  and  a 
Hamiltonian  function  H(p\, . . .  ,p„;  q\, . . . ,  qn\  i).  By  Darboux's  theorem,  one 
can  choose  a  local  coordinate  system  (pi, . . .  ,pn;  q\, . . . ,  qn)  such  that  the  sym- 
plectic two-form  has  the  form 

n 
u2  =  ^T  dpi  A  dqi.  (2.1) 

i=l 

Then,  by  Stoke's  lemma,  a  finite-dimensional,  unconstrained  Hamiltonian  sys- 
tem can  by  characterized  by  Hamilton's  equations  of  motion: 

dqi_=dH_        dpi_      _dB_ 

dt        dp,'        dt  dqi'  l      ' 

where  p,  and  <j,  are  the  phase  space  coordinates,  position  and  momentum. 

Consider  a  Hamiltonian  dynamical  system  with  N  degrees  of  freedom;  this 
system  has  a  2N-dimensional  phase  space.  A  Hamiltonian  system  is  integrable 
if  it  has  at  least  N  independent  constants  of  the  motion,  including  the  energy, 
for  which,  for  any  two  such  constants  F\  and  F2,  the  Poisson  bracket 

dF1dF2      dF2dF1 

{Fl>F2]--dq-dp---dq-dp-  <L13) 

must  vanish.  The  constants  F\  and  F2  are  then  said  to  be  in  involution. 
Any  function  which  is  in  involution  with  the  Hamiltonian  is  a  constant  of  the 
motion,  which  can  be  seen  using  Hamilton's  equations. 
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For  a  time-independent,  one  degree-of-freedom  Hamiltonian  system,  since 
the  energy  is  a  constant  of  the  motion,  we  have  H(p,  q)  =  E,  so,  up  to  sign 
ambiguities,  the  momentum  can  be  written  as  a  function  of  the  position  alone, 
P  —  Pi<l,E)-  With  the  help  of  Hamilton's  equations,  we  can  then  reduce  to 
quadratures 


t 


-I 


dq 


(dH/dp) 


(2.3) 
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Thus,  the  system  is  integrable,  though  we  may  not  be  able  to  evaluate  the 
integral  analytically. 

Trajectory 


Figure     2.1     Action-angle  variables.    The  figure  shows  a  trajectory  on  the 
surface  of  a  torus  for  a  two  degree-of-freedom  integrable  Hamiltonian. 


If  a  Hamiltonian  system  is  integrable,  action-angle  variables  can  be  found. 
The  action  variables  J,  are  defined  to  be 


Ji 


i-^fpidqi, 


(2.4) 
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and  the  angles  0,  are  the  conjugate  coordinates  to  Jt.  In  this  case,  motion  lies 
on  an  N-dimensional  torus,  for  which  the  radii  of  the  torus  are  determined  by 
the  actions.  As  a  consequence  of  this  transformation,  the  Hamiltonian  can  be 
written  as  a  function  of  the  actions  only,  i.e.  H  =  H(3);  the  angles  are  cyclic 
coordinates.  The  equations  of  motion  are  then  seen  to  be  trivial: 

Mi      dH  dJi         dH 

it-m"*'  -di=-wr^  (25) 

where  the  Ui  are  the  frequencies,  so  that  the  complete  solution  to  the  equations 
of  motion  is  simply 

Bi  =  uit  +  (j>i.  (2.6) 

For  such  systems  the  coordinates  p,  and  <?,  can  be  written  as  a  Fourier  expansion 
of  the  angle  variables.  This  is  said  to  yield  multiperiodic  (or  quasiperiodic) 
motion  in  pi  and  ft. 

Now  consider  a  two  degree-of-freedom  integrable  system,  for  which  the  mo- 
tion lies  on  a  two-dimensional  torus  (see  figure  2.1).  If  the  two  frequencies  are 
commensurate  (rational  multiples  of  one  another),  then  the  motion  is  periodic 
and  the  orbit  will  repeat  itself.  However,  in  the  generic  case,  the  frequencies  are 
irrational  multiples  of  one  another,  and  the  orbit  will  cover  the  entire  surface 
of  the  torus.  These  tori  in  action-angle  space  are  called  invariant  tori. 

2.2  Near-integrable  Systems 
Near-integrable  systems  are  often  described  as  perturbations  of  integrable 
systems.   The  Hamiltonian  for  such  a  system  can  be  written  in  terms  of  the 
action-angle  variables  of  the  unperturbed  motion, 

H(J,6)  =  Ho(J)  +  eH1(J,0),  (2.7) 
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where  Ho  is  the  unperturbed,  integrable  Hamiltonian,  and  e  is  a  small  parame- 
ter. Even  though  integrability  is  an  exceptional  circumstance,  many  dynamical 
systems  are  known  that  are  'close  to  integrable,'  for  example,  the  three-body 
problem  and  the  fast  rotations  of  a  rigid  body. 

In  order  to  solve  equation  (2.7),  one  has  to  resort  to  perturbation  theory. 
The  conventional  approach  is  to  expand  the  Hamiltonian  in  a  Fourier  series  in 
the  parameter  e.  A  generating  function  can  be  introduced  to  transform  to  a 
new  set  of  action-angle  variables.  The  Hamiltonian  for  this  new  set  of  variables 
can  be  expanded  in  a  Taylor  series  in  e.  The  consequence  of  carrying  out  this 
procedure  is  that  the  actions  are  now  expressed,  to  lowest  order  in  e,  with  a 
denominator  which  may  be  small.  This  denominator  can  be  expressed  as  a 
term  m  •  w,  which  leads  to  resonances  when  the  values  of  m  are  such  that  the 
denominator  approaches  zero.  As  was  originally  shown  by  Poincare  [16],  such 
perturbation  expansions  are  expected  to  diverge  near  these  internal  resonances. 

Kolmogorov  [17],  in  an  attempt  to  overcome  the  problem  of  small  denom- 
inators associated  with  resonance  zones,  constructed  a  superconvergent  per- 
turbation theory  which  was  applicable  to  nonresonant  tori  .  For  integrable 
systems,  the  phase  space  is  completely  divided  into  invariant  tori.  For  a  sys- 
tem that  is  close  to  an  integrable  system,  Kolmogorov  proved  that  the  motion 
remains  quasi-periodic  for  most  initial  data.  Nonresonant  tori,  or  KAM  tori 
named  after  Kolmogorov,  Arnol'd  and  Moser,  are  not  destroyed  by  resonances 
when  a  sufficiently  small  parameter  e  is  introduced.  Instead,  phase  space  is 
divided  into  regions  where  invariant  tori  exist,  and  regions  where  resonances 
have  destroyed  the  tori  and  motion  is  stochastic.  The  volume  of  phase  space 
occupied  by  resonances  approaches  zero  as  e  approaches  zero.  As  e  increases, 
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Figure  2.2  The  figure  shows  the  intersection  of  a  periodic  orbit  with  a 
Poincare  surface  of  section.  The  resulting  Poincare  map  consists  of  a  period- 
three  orbit. 


the  measure  of  the  stochastic  region(s)  increases,  as  more  and  more  KAM  tori 
are  destroyed.  These  ideas  are  more  clearly  seen  after  we  introduce  the  concept 
of  Poincare  surfaces  of  section,  which  is  the  subject  of  the  next  section. 


2.3  Poincare  Maps 
Poincare  had  the  useful  idea  to  convert  the  study  of  continuous  systems  to 
that  of  discrete  systems  by  reducing  the  number  of  dimensions  of  the  system. 
The  resulting  system  is  called  a  Poincare  map,  or  surface  of  section.  Consider 
an  autonomous  Hamiltonian  system  with  two  degrees-of- freedom,  i.e.,  with 
a  four-dimensional  phase  space.  Any  trajectory  lies  on  a  three-dimensional 
energy  hypersurface  H(pi,P2,qi,qi)  =  #0,  so  p2  =  P2(pi, 91,92).  We  can  then 
study  the  crossings  of  the  trajectory  across  a  two-dimensional  surface  given, 
say,  by  qi  =  const.,  in  a  certain  direction,  say,  q\  >  0. 
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If  the  system  is  integrable,  another  global  constant  of  the  motion  exists, 
I(pi,P2,<li,<l2))  —  const.,  and  so  we  can  write  p\  =  Pi(qi,q2)-  In  this  case 
successive  crossing  lie  on  regular  curves  on  the  surface  of  section. 

If  a  second  integral  does  not  exist,  intersections  of  the  trajectory  with  the 
surface  of  section  do  not  necessarily  lie  on  regular  curves.  Consider,  as  an 
example,  the  Henon-Heiles  potential: 

U{x,y)=l-{x2  +  y2  +  2x2y-l-y*).  (2.8) 

Figure  2.3  is  a  surface  of  section  for  this  potential  computed  for  three  different 
values  of  the  energy  [18].  For  the  case  E  =  1/12,  it  appears  that  all  trajectories 
lie  on  smooth  curves.  Therefore,  this  potential  may  appear  to  be  integrable 
for  this  value  of  the  energy,  even  though  it  is  not.  As  one  increases  the  energy, 
however,  these  regular  curves  appear  to  break  down.  One  can  see  that,  for 
E  =  1/8,  the  surface  is  partially  covered  by  invariant  tori,  and  part  of  the 
phase  space  is  stochastic.  For  E  =  1/6,  most  of  the  phase  space  is  stochastic. 
The  construction  of  a  Poincare  map  yields  several  advantages  for  the  study 
of  differential  equations,  as  noted  by  Wiggins  [19].  (1)  Dimensional  reduction 
-  The  construction  of  a  Poincare  map  removes  at  least  one  of  the  variables  of 
the  original  equation.  This  also  means  that  some  of  the  information  about  the 
original  system  is  lost.  (2)  Global  dynamics  -  For  systems  with  few  degrees-of- 
freedom,  computation  of  a  Poincare  map  may  provide  insight  into  the  global 
dynamics  of  the  system.  Unfortunately,  the  choice  of  cross-section  may  be  of 
importance;  different  surfaces  of  section  for  the  same  system  may  not  necessar- 
ily have  the  same  dynamics.  However,  there  are  situations  when  the  Poincare 
map  contains  all  the  information  regarding  the  dynamics  of  the  original  dif- 
ferential equation,  as,  for  example,  when  the  region  of  phase  space  of  interest 
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Figure    2.3     Poincare  surface  of  section  for  the  Henon-Heiles  potential 


(a)  E  =  l/12,(b)  E  =  1/8,  (c)  E  =  1/6 
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contains  only  periodic  orbits.  (3)  Conceptual  clarity  -  A  fixed  point  in  the 
surface  of  section  corresponds  to  a  periodic  orbit  of  the  original  differential 
equation.  It  can  be  shown  that  if  a  fixed  point  of  the  Poincare  map  is  asymp- 
totically stable,  then  the  corresponding  periodic  orbit  is  also  asymptotically 
stable.  It  is  easier  to  study  the  stability  of  the  fixed  point,  which  is  simply 
characterized  in  terms  of  the  eigenvalues  of  the  map  linearized  about  the  fixed 
point,  than  the  stability  of  the  periodic  orbit  in  the  original  differential  equa- 
tion. Results  concerning  the  dynamics  of  the  differential  equation  may  be  more 
clearly  stated  for  the  Poincare  map. 

Moreover,  it  is  apparent  that  maps  are  much  easier  to  study  numerically 
than  differential  equations.  Two-dimensional  maps  that  are  dynamically  simi- 
lar to  the  surfaces  of  section  of  nonintegrable  Hamiltonians  can  be  easily  con- 
structed, as,  e.g.,  the  figures  associated  with  the  Henon-Heiles  potential.  Thus, 
studying  the  dynamics  of  such  maps  should  yield  useful  information  regarding 
the  dynamics  of  differential  equations.  For  simplicity,  it  would  be  preferable  to 
learn  as  much  as  possible  about  two-dimensional  Hamiltonian  mappings  before 
tackling  the  more  difficult  task  of  studying  Hamiltonian  differential  equations, 
with  the  understanding  that  such  a  mapping  should  be  qualitatively  similar  to 
a  Poincare  map. 

2.4  Invariant  Manifolds  and  Homoclinic  Points 
Invariant  manifolds,  including  stable  and  unstable  manifolds,  are  important 
analytical  tools  for  the  study  of  dynamical  systems.  A  set  S  is  said  to  be 
invariant  if,  under  the  action  of  the  vector  field  or  map,  a  point  within  S 
remains  within  S  for  all  time.  As  far  as  we  are  concerned,  a  manifold  can  be 
considered  to  be  a  set  which  locally  has  the  structure  of  Euclidean  space.  In 
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Figure  2.4  The  figure  shows  one- dimensional  stable  and  unstable  manifolds 
of  a  fixed  point,  7.  The  stable  and  unstable  manifolds  of  this  point,  Ws  and 
Wu,  are  tangent  to  the  linear  subspaces  Es  and  Eu. 


nonlinear  settings,  it  is  sufficient  to  consider  a  manifold  to  be  an  m-dimensional 
surface  embedded  in  Rn. 

Note  that  when  we  call  something  a  stable  or  unstable  manifold,  it  has 
to  be  of  something.  Most  often,  these  manifolds  are  considered  to  be  of  fixed 
points.  However,  invariant  manifolds  of  more  general  sets,  e.g.,  typical  stochas- 
tic trajectories  in  Hamiltonian  systems,  may  also  be  of  importance. 

For  fixed  points,  the  stable,  unstable  and  center  manifolds  span  the  sub- 
spaces  determined  by  their  eigenvalues,  sufficiently  close  to  that  fixed  point; 
positive  eigenvalues  correspond  to  unstable  manifolds,  negative  to  stable  man- 
ifolds, and  zero  to  center  manifolds.  For  fixed  points,  the  matrix  of  first  partial 
derivatives  is  constant.  For  more  general  points,  however,  this  is  not  true  and 
calculation  of  the  manifolds  are  better  done  numerically  (by  linearizing,  as  for 
fixed  points).  Note  that,  as  illustrated  in  figure  2.4,  the  nonlinear  manifolds  are 
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Figure  2.5  The  figure  shows  one-  and  two-dimensional  stable  and  unstable 
linear  subspaces  of  a  fixed  point.  Nearby  trajectories  tend  to  diverge  from  the 
fixed  point  in  the  unstable  directions  and  converge  towards  the  fixed  point  in 
the  stable  directions. 


tangent  to  the  associated  linear  manifolds  at  the  location  of  the  fixed  point. 
Figure  2.5  [20]  shows  several  one-  and  two-dimensional  stable  and  unstable 
subspaces  and  the  behavior  of  nearby  trajectories. 

Consider  the  stable  and  unstable  manifolds  of  a  hyperbolic  periodic  point 
7.  A  hyperbolic  point  is  a  point  for  which  the  manifolds  intersect  transversely, 
that  is,  non-tangentially.  It  is  sufficient  to  consider  the  periodic  point  to  be 
a  fixed  point,  simply  by  applying  any  argument  for  /  to  fp,  where  /  is  the 
given  map,  and  fp  is  the  pth  iterate  of  the  map.  A  point  that  is  within  the 
intersection  of  the  manifolds  is  said  to  be  homoclinic  to  7.    The  existence  of 
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Figure     2.6     The  figure  shows  a  homoclinic  tangle  caused  by  the  infinite 
number  of  intersections  of  the  stable  and  unstable  manifolds  of  a  fixed  point. 


such  points  is  an  important  'trigger'  for  chaotic  behavior.  It  can  also  be  shown 
that  the  existence  of  a  single  such  point  implies  infinitely  many  such  points. 
This  is  called  the  homoclinic  tangle,  and  is  shown  in  figure  2.6.  Moreover, 
transverse  intersections  of  the  stable  and  unstable  manifolds  are  a  necessary 
condition  for  structural  stability.  Difficulty  occurs  when  the  intersections  are 
not  transverse.  When  they  are  tangent,  we  obtain  a  homoclinic  tangency,  and 
the  assumption  of  hyperbolicity  is  lost.  Regions  where  such  a  tangency  occurs 
are  then  nonhyperbolic.  Since  physical  systems  are,  in  general,  nonhyperbolic, 
the  behavior  of  the  stable  and  unstable  manifolds  is  of  vital  importance  in 
understanding  the  structural  stability  of  dynamical  systems. 


2.5  Liapunov  Exponents  and  the  Definition  of  Chaos 
Liapunov  exponents  provide  an  important  tool  in  the  understanding  of  the 
stochastic  regions  of  phase  space.  Nearby  orbits  in  stochastic  regions  diverge 
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from  each  other  exponentially.    This  is  called  sensitive  dependence  on  initial 
conditions,  and  is  the  fundamental  phenomenon  associated  with  chaotic  behav- 
ior. Liapunov  characteristic  exponents  are  a  primary  quantity  in  determining 
the  existence  and  rate  of  this  divergence. 

Consider  two  trajectories  with  initial  conditions  x(0)  and  x(0)  +  £x(0), 
respectively.  Exponential  sensitivity  to  initial  conditions  would  imply  that,  for 
a  long  enough  time  t,  6x(t)  ~  exp(A£)£x(0),  for  some  positive  A.  In  a  regular 
region  where  exponential  sensitivity  does  not  occur,  one  would  expect  that 
A  — ►  0  as  t  — *•  oo. 

The  Liapunov  exponents  can  be  explicity  defined  as  follows  [21,22].  Given 
an  initial  condition  x(0)  and  an  initial  seperation  vector  £x(0),  we  can  define 
the  mean  exponential  rate  of  divergence  of  the  trajectories  to  be 

where  d(x(0),  t)  is  simply  the  norm  of  6x.(t),  usually  chosen  to  be  the  Euclidean 
norm,  though  it  is  of  course  independent  of  the  choice  of  norm. 
One  can  then  find  an  orthogonal  basis  {e,}  of  6x  such  that 

A,(x(0))  =  A(x(0),e,),  (2.10) 

where  the  Aj  are  the  Lyapunov  exponents.  Once  the  largest  Liapunov  exponent 
is  found,  the  others  can  be  found  using  an  orthogonalization  procedure,  such  as 
that  of  Gram-Schmidt.  The  largest  Liapunov  exponents  for  numerous  initial 
conditions  in  the  Henon-Heiles  potential  were  calculated  by  Benettin  et  al. 
[21].  Several  of  these  exponents  are  re-calculated  here,  and  shown  in  figure  2.7. 
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Figure  2.7  The  figure  shows  the  largest  Liapunov  exponent  for  six  sets  of 
initial  conditions  in  the  Henon-Heiles  potential,  with  E  =  0.125.  Three  of 
these  are  within  the  stochastic  region  and  seem  to  converge  to  the  same  value, 
while  the  other  three  are  within  the  regular  region  and  tend  towards  zero. 


For  Hamiltonian  systems,  Liouville's  theorem  implies  that  the  number  of 
positive  Liapunov  exponents  equals  the  number  of  negative  ones.  Hence  two- 
dimensional  autonomous  Hamiltonian  systems  have  only  one  positive  Liapunov 
exponent. 

Despite  the  fact  that  an  exact  definition  of  chaos  is  not  universally  agreed 
upon,  the  existence  of  positive  Liapunov  exponents  is  critical  to  any  reasonable 
definition.  Devaney's  topological  definition  of  chaos  is  frequently  cited  [23]. 
According  to  this  definition,  there  are  three  properties  all  chaotic  systems 
must  possess:  (1)  a  dense  set  of  periodic  orbits;  (2)  topological  transitivity, 
the  property  that  an  orbit  beginning  in  one  arbitrarily  small  neighborhood 
will  eventually  pass  through  any  other;  and  (3)  sensitive  dependence  on  initial 
conditions.    It  has  actually  recently  been  shown  that,  generically,  conditions 
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(1)  and  (2)  imply  condition  (3)  [24]  .  It  would  be  possible,  however,  to  find 
a  system  which  satisfies  condition  (3)  while  not  obeying  (1)  and  (2),  although 
such  a  system  would  be  pathological. 

Reichl  [25]  defines  a  chaotic  system  to  be  a  system  with  positive  Kol- 
mogorov  entropy,  the  sum  of  the  positive  Liapunov  exponents.  Such  systems 
are  called  K-systems.  A  more  restricted  set  of  systems  are  C-systems,which 
are  defined  to  be  systems  that  have  a  decomposable  linear  tangent  space,  i.e., 
these  are  systems  that  are  everywhere  hyperbolic.  Such  systems  are  not  very 
physical  in  the  sense  that  this  suggests  that  the  Hamiltonian  is  discontinuous. 
These  hyperbolic  systems  may  be  considered  the  most  chaotic,  in  the  sense 
that,  for  any  choice  of  initial  conditions,  one  has  the  maximum  exponential 
sensitivity.  In  any  case,  calculating  a  positive  Liapunov  exponent  is  a  good 
indication  that  chaos  is  present  in  a  given  system. 

2.6  Dissipative  Systems 
Most  often,  chaos  is  considered  in  the  realm  of  dissipative  systems,  which 
are  characterized  by  the  presence  of  attractors  in  the  phase  space.  Chaotic  dis- 
sipative systems  have  strange  attractors,  on  which  trajectories  have  sensitive 
dependence  on  initial  conditions.  Conservative  dynamical  systems,  i.e.,  Hamil- 
tonian systems,  have  no  such  attractors  because  trajectories  in  the  phase  space 
are  volume  preserving  -  attractors  are  characterized  by  a  reduction  of  dimen- 
sion. Nevertheless,  regions  of  the  phase  space  in  nonintegrable  Hamiltonian 
systems  are  often  'chaotic'  in  this  sense  that  they  have  sensitive  dependence 
on  initial  conditions.  Such  systems  are  often  called  stochastic,  as  their  phase 
space  is  divided  into  'regular'  and  'stochastic'  regions. 
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Dissipative  dynamical  systems  can  be  studied  in  much  the  same  way  as 
Hamiltonian  systems  are.  The  standard  methods  of  determining  Poincare  sur- 
faces, stable  and  unstable  manifolds,  and  Liapunov  exponents  still  apply.  Many 
systems  of  physical  interest  are  best  modelled  by  dissipative  dynamical  systems. 
However,  nature  may  be  Hamiltonian  at  its  most  fundamental  level,  and  the 
complex  dynamics  that  arise  in  nonintegrable  Hamiltonian  systems  make  these 
worthy  of  study  in  themselves. 


CHAPTER  3 

CHAOS  IN  THE  GRAVITATIONAL 

N-BODY  PROBLEM 

In  general,  the  problem  of  N  bodies  interacting  via  Newton's  law  of  gravita- 
tion is  nonintegrable  for  TV  >  2.  The  full  6iV-dimensional  problem,  in  fact,  has 
ten  integrals  of  motion,  related  to  the  initial  conditions  of  the  center  of  mass 
(3),  the  total  linear  momentum  (3),  the  total  angular  momentum  (3),  and  the 
energy  (1).  The  full  JV-body  system  has  been  shown  numerically  to  be  chaotic 
in  the  sense  that  it  has  positive  Lyapunov  exponents.  Consider  evolving  an 
iV-body  system  in  its  6./V-dimensional  phase  space,  and  evolving  a  second  such 
system  which  is  very  nearly  identical,  i.e.  the  system's  location  in  phase  space 
is  very  close  to  that  of  the  first  system.  Numerical  studies  indicate  that,  for 
N  >  2,  these  two  trajectories  will  diverge  from  each  other  exponentially,  which 
is  an  indication  of  the  presence  of  chaos. 

We  are  of  course  primarily  interested  in  the  case  when  N  >  2,  for  the 
case  of  galactic  dynamics.  For  the  study  of  such  systems,  numerous  numerical 
codes  exist,  a  standard  one  being  Aarseth's,  found,  e.g.,  in  Appendix  4.B  of 
Galactic  Dynamics  [2].  The  primary  purpose  of  such  codes  is  to  integrate  the 
equations  of  motion  directly, 

..        -dU 
mi9i  =  — ,  (3.1) 

where 

JT  _  ^-\         Gmirrij 

^k^H^'  (3'2) 
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Computationally,  equation  (3.2)  leads  to  problems,  since  the  denominator 
blows  up  when  q,  is  close  to  q,.  The  standard  method  of  dealing  with  this 
problem  is  to  insert  a  'softening  parameter'  into  the  equation,  so  that  the 
softened  potential  reads 

°~  E  £&.  (") 

i<«£i<ivv     + 
where  e  is  the  softening  parameter,  and  r  =||  q,-  —  qy  ||. 

The  softening  parameter  can  also  be  interpreted  as  'turning  off'  the  effects 
of  close  encounters  of  nearby  particles  in  the  TV-body  simulation.  One  can  thus 
determine  the  effect  that  this  close  encounter  'noise'  has  on  the  instability  of 
the  gravitational  TV-body  problem,  compared  to  the  effect  of  the  bulk  mean 
field  of  the  system.  In  this  chapter,  this  effect  will  be  studied  by  performing 
TV-body  simulations  to  determine  the  time-scale  of  instability  for  various  values 
of  the  softening  parameter. 

3.1  Method 

This  section  describes  the  method  by  which  the  time-scale  of  instability  for 
each  simulation  is  determined.  The  simulations  were  run  using  an  TV-body  code 
developed  by  Haywood  Smith,  which  is  similar  to  the  basic  code  of  Aarseth. 

For  each  case,  the  particles  were  distributed  randomly  with  uniform  prob- 
ability within  a  sphere  of  radius  R0  =  1,  with  velocities  chosen  similarly,  but 
scaled  such  that  the  initial  virial  ratio  of  the  system,  -2Et/Ep  =  1.  The 
particles  then  interacted  via  forces  derived  from  equation  (3.3).  A  perturbed 
system  was  then  generated  in  an  identical  fashion,  except  that  each  particle 
was  displaced  by  a  distance  Sr  =  2  x  10~8  from  the  unperturbed  particles  in  a 
randomly  chosen  direction. 
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These  simulations  were  run  for  systems  comprised  of  N  =  360  particles,  300 
of  which  had  a  mass  m  =  1  (light  particles)  and  the  remaining  60  of  which  had 
m  =  5  (heavy  particles).  In  addition,  some  simulations  were  run  for  systems 
comprising  of  N  =  996  particles,  with  the  same  ratio  of  light/heavy  particles. 
Three  different  runs  were  performed  for  each  value  of  c,  and  the  results  from 
these  runs  were  averaged  to  extract  mean  properties. 

For  each  pair  (perturbed  and  unperturbed)  of  simulations,  the  perturbation 
of  each  particle's  position,  <5r,(<),  was  determined,  as  well  as  the  perturbation 
in  velocity,  6v{(t).  These  data  were  then  analyzed  to  determine  two  types  of 
quantities. 

First,  the  time  taken  for  each  perturbation  6r,(r)  to  grow  to  a  value  e  times 
the  initial  6r,(0)  =  2  x  10~8  was  determined.  These  times  were  then  analyzed 
to  extract  mean  e-folding  times  t*ti  and  r+)#  for  the  light  and  heavy  particles, 
as  well  as  the  standard  deviations  o+i  and  0+ jj.  To  ensure  that  the  qualitative 
results  were  not  dependent  on  the  choice  of  e  as  an  amplification  factor,  an 
identical  analysis  was  also  effected  for  the  1.2-folding  times. 

The  total  3JV-dimensional  configuration  and  velocity  space  perturbations 
were  computed  as  well.  These  are  defined  to  be 

Ar*|b(J>rip)      and      Av  =  iln(£  |*v,|2).  (3.4) 

i  » 

The  corresponding  quantities  Ar^,  Avi,  Arjj  and  Avjj,  associated  individually 
with  the  light  and  heavy  particles,  were  also  computed.  Growth  time  scales 
iA,i  and  <a  h  were  then  extracted,  derived  from  the  quantities  Ari  and  Arjj 


via 


Ar(0  =  Ar(0)  +  (</*A).  (3.5) 
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Characteristic  growth  times  for  the  velocity  perturbations  were  determined  as 
well,  which  agreed  with  the  time  scales  extracted  from  the  configuration  space 
perturbations  but,  owing  to  the  larger  irregularities,  exhibited  larger  statistical 
uncertainties. 

3.2  Results 

Table  3.1 

Average  e-folding  times  and  Standard  Deviations  (N  =  360) 


€ 

t*,LftcT 

i*,Hftcr 

&*,Lltcr 

°*,Hftcr 

lO"8 

0.47  ±0.01 

0.49  ±  0.04 

0.21  ±0.01 

0.21  ±  0.01 

0.005 

0.50  ±  0.01 

0.50  ±  0.05 

0.23  ±0.01 

0.21  ±0.01 

0.01 

0.54  ±0.01 

0.55  ±  0.03 

0.23  ±0.01 

0.22  ±0.01 

0.02 

0.66  ±  0.02 

0.66  ±  0.06 

0.26  ±  0.01 

0.24  ±0.02 

0.05 

1.09  ±0.02 

1.12  ±0.05 

0.41  ±  0.02 

0.41  ±  0.05 

0.08 

1.71  ±0.03 

1.79  ±0.08 

0.56  ±  0.03 

0.54  ±  0.02 

0.08 

2.65  ±  0.03 

2.49  ±  0.35 

0.90  ±0.02 

0.95  ±0.07 

0.08 

3.58  ±  0.25 

3.79  ±  0.30 

1.02  ±0.16 

1.00  ±  0.22 

The  mean  e-folding  times  for  both  light  and  heavy  particles,  and  the  asso- 
ciated dispersions,  are  tabulated  in  Table  3.1  for  the  different  values  of  e.  It 
can  be  seen  immediately  that  the  e-folding  times  t+yL  and  t^iH  for  fixed  e  are 
in  fact  equal  to  within  statistical  uncertainties.  In  terms  of  the  first  two  mo- 
ments, there  are  no  discernable  differences  between  light  and  heavy  particles. 
This  similarity  extends  also  to  the  full  distributions  of  e-folding  times,  which 
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Figure     3.1     The  mean  e- folding  times  for  light  and  heavy  particles  as  a 
function  of  e/l,  the  ratio  of  the  softening  parameter  to  the  interparticle  spacing. 
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were  binned  and  examined  graphically.  There  is  no  evidence  whatsoever  that 
the  distributions  for  the  light  and  heavy  particles  are  not  identical. 

It  is  also  observed  that,  for  e  <  0.18,  the  values  of  the  time  scales  t,  i 
and  t^n  are  smooth,  monotonically  increasing  functions  of  e.  In  fact,  these 
values  are  consistent  with  a  power-law  growth  with  respect  to  the  soften- 
ing parameter,  if  it  is  expressed  in  units  of  the  typical  interparticle  spacing, 
I  =  T(4/3)iV-1/3jR0  m  0.890JV-V3,  computed  following,  e.g.,  Chandrasekhar 
(1941).  These  data  are  easily  fit  to 


u(€/i)  =  u(o)  +  M*W, 


(3.6) 
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Figure     3.2     The  mean  e- folding  times  for  light  and  heavy  particles  as 
function  of  e/l,  plotted  on  a  linear  scale. 


1.5 


with  best  fit  values  pL  =  1.35  ±  0.08  and  Ai  =  2.15  ±  0.03<cr  for  the  light 
particles  and  pjj  =  1.55  ±  0.23  and  Ah  =  2.35  ±  0.16tCT  for  the  heavy  particles. 
If  the  point  e  =  0.18  be  excluded  from  the  fit,  the  values  of  the  exponents 
are  increased  to  pi  =  1.38  and  pH  =  1.60.  The  goodness  of  fit  is  illustrated 
in  figure  3.1.  Because  there  are  only  1/5  as  many  heavy  particles  as  light 
particles,  the  uncertainties  in  t+tH  are  quite  large,  so  large  that  p#  and  pi  are 
in  agreement  within  statistical  uncertainties.  Indeed,  if  one  excludes  the  value 
e  =  0.005,  the  best  fit  value  of  pjj  reduces  to  1.37. 
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Figure    3.3     The  mean  1.2-folding  times  for  light  and  heavy  particles  as  a 
function  of  e/l. 


The  values  for  the  e-folding  times  can  also  be  fit  by  a  simple  linear  rela- 
tionship if  values  e  <  0.02  are  excluded.  The  best  fit  values  for 

U{t/l)  =  U(0)  +  £♦(«//)  (3.7) 

are  then  <*,.t(0)  =  0.24*cr  and  B^L  =  2.38*cr  for  the  light  particles  and 
t*,H(Q)  =  0.19tcr  and  £?,,#  =  2.48<cr  for  the  heavy  particles.  The  goodness 
of  fit  is  illustrated  in  figure  3.2. 

The  1.2-folding  times  th2,L  and  h.2,H ,  plotted  in  Figure  3.3,  are  also  well  fit 
by  a  simple  power  law  growth.  In  this  case,the  best  fits  are  p\.2,L  =  1.50±0.12 
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and  A\.2,l  =  1.54±0.05tfcr  for  the  light  particles  and  p\.2,H  =  1.50±0.18  and 
•^1.2, H  =  1.54±0.08tfcr  for  the  heavies.  The  data  are  thus  consistent  with  the 
interpretation  that  the  ratios  ti.2,L/i*,L  and  ti.2,B/t*,B  aie  independent  of  e, 
although  there  it  a  hint  that  the  1.2-folding  times  may  grow  somewhat  more 
rapidly.  For  both  the  light  and  heavy  particles,  the  best  fit  value  for  the  ratio 
is  <i.2/**  =  0.61  ±  0.03. 

Table  3.2 

Growth  Times  (JV  =  360) 


e 

*A,l/*cr 

*A,ir/*er 

io-8 

0.10  ±0.01 

0.09  ±0.01 

0.005 

0.15  ±0.01 

0.16  ±0.01 

0.01 

0.18  ±0.02 

0.18  ±0.01 

0.02 

0.28  ±0.02 

0.27  ±  0.03 

0.05 

0.52  ±0.02 

0.52  ±  0.07 

0.08 

0.80  ±  0.03 

0.82  ±  0.04 

0.08 

1.13  ±0.06 

1.27  ±0.03 

0.08 

1.17  ±0.21 

1.34  ±0.25 

The  mean  growth  times  t\  for  the  total  perturbations  Ari  and  Arjj  of 
the  light  and  heavy  particles  are  tabulated  in  Table  3.2.  Again,  no  difference 
can  be  seen  between  the  results  for  the  light  and  heavy  particles,  to  within 
statistical  uncertainties. 

These  growth  times  are  smooth,  monotonically  increasing  functions  of  e  as 
well.  As  illustrated  in  figure  3.4,  they  are  well  fit  by  a  power  law 


tA(t/l)  =  tA(0)  +  AA(e/iy. 


(3.8) 
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Figure     3.4     The  growth  times  for  the  total  perturbations  as  a  function  of 
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The  best  fit  values  here  are  j>a,L  =  1.13  ±  0.05  and  A\tL  =  1.12  ±  0.04<cr  for 
the  light  particles  and  pA>Jsr  =  1-12  ±  0.06  and  Aa,h  =  1-23  ±  0.05*cr  for  the 
heavy  particles. 

Several  simulations  were  run  for  the  case  N  =  996  as  well.  We  had  assumed 
previously  that  the  instability  time  scales  depended  only  on  the  value  of  the 
ratio  e/l.  The  mean  interparticle  spacings  for  the  systems  with  N  =  360  and 
N  =  996  differ  by  a  factor  of  (360/996)1/3  «  0.712.    If  our  assumption  is 
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Table  3.3 

Average  e-folding  times  and  Standard  Deviations  (TV  = 


360) 


e 

N 

**,LI*cr 

U,u/icr 

CT*,L/tcr 

<T*,H/tcr 

lO"8 

360 

0.47  ±  0.01 

0.49  ±  0.04 

0.21  ±  0.01 

0.21  ±0.01 

lO"8 

996 

0.45  ±  0.01 

0.46  ±  0.02 

0.19  ±0.01 

0.19  ±0.01 

0.0200 

360 

0.66  ±  0.02 

0.66  ±  0.06 

0.26  ±  0.01 

0.24  ±  0.02 

0.0143 

996 

0.67  ±0.01 

0.67  ±0.03 

0.26  ±  0.02 

0.25  ±  0.02 

h.2,Lh<* 

tl.2,H/tcr 

0-*,iAcr 

<7"*,lAcr 

*A,lAcr 

t\,H/tcr 

0.31  ±0.01 

0.31  ±  0.03 

0.19  ±0.01 

0.18  ±0.01 

0.10  ±0.01 

0.09  ±0.01 

0.30  ±0.01 

0.30  ±  0.01 

0.18  ±0.01 

0.18  ±0.01 

0.12  ±0.02 

0.12  ±0.02 

0.40  ±  0.01 

0.38  ±  0.03 

0.24  ±  0.01 

0.22  ±  0.02 

0.28  ±  0.03 

0.27  ±0.05 

0.40  ±  0.01 

0.41  ±0.01 

0.23  ±0.01 

0.22  ±  0.01 

0.28  ±  0.04 

0.30  ±  0.03 

correct,  the  time  scales  should  be  identical  for  these  two  systems  by  decreasing 
the  softening  parameter  by  a  factor  of  0.712  for  the  case  of  996  particles. 

The  results  for  996  particles  are  tabulated  and  compared  to  the  results 
for  360  particles  in  table  3.3,  for  two  values  of  the  softening  parameter.  In 
these  cases,  the  predicted  scaling  is  confirmed,  and  it  can  also  be  seen  visually 
that  the  shapes  of  the  distributions  are  identical.  Thus,  the  effect  of  softening 
seems  to  be  independent  of  the  number  of  particles,  which  implies  that  close 
encounters  remain  important  even  for  large  N. 


3.3  Representative  Particle  Trajectories 
The  form  of  the  individual  particle  perturbations  \6r\  and  \8v\  and,  partic- 
ularly, the  question  of  the  degree  to  which  different  particles  behave  similarly, 
was  studied  by  extracting  and  then  examining  individually  the  perturbations 
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Figure   3.5     The  seperation  of  typical  trajectories  |<$r(tf)|  for  three  individual 
particles  each  with  c  =  0.005,0.05,  and  0.18. 


of  5%  each  of  the  light  and  heavy  particles  in  several  different  simulations.  It 
was  immediately  evident  that  there  are  no  significant  differences  between  the 
trajectories  of  the  light  and  heavy  particles,  and,  for  this  reason,  in  the  fol- 
lowing no  attempt  has  been  made  to  distinguish  between  the  light  and  heavy 
particles.  This  analysis  of  individual  particles  was  effected  for  three  different 
values  of  epsilon,  c  =  0.005,  0.05,  and  0.18.  Three  runs  each  were  analyzed  for 
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e  =  0.005  and  0.18,  but  only  one  run  was  analyzed  for  the  intermediate  value 
c  =  0.05.  The  behavior  of  three  representative  particles  for  each  of  these  values 
of  e  is  illustrated  in  Figure  3.5. 

It  is  evident  that,  regardless  of  the  value  of  e,  there  is  a  roughly  exponential 
growth  for  the  individual  particle  perturbations,  interwoven,  however,  with 
significant  "bumps."  This  overall  qualitative  behavior  can  be  interpreted  as 
representing  a  sum  of  two  effects,  namely  (1)  a  smooth  growth,  which  can 
be  attributed  to  the  bulk  mean  field,  and  (2)  sudden,  random  kicks  which 
can  be  attributed  to  shorter  range  interactions  between  particularly  proximate 
particles.  Significantly,  however,  there  is  an  initial  transient  period,  lasting 
about  0.5tfcr  to  licr,  during  which  the  perturbation  does  not  exhibit  an  overall 
rapid  exponential  growth.  In  fact,  the  perturbation  tends  initially  to  either 
decrease  somewhat  or  to  only  grow  relatively  slowly.  A  more  rapid  growth 
does  not  appear  to  set  in  until  it  is  triggered  by  some  initial  shorter  range 
encounter. 

For  all  values  of  c,  the  individual  velocity  perturbations  \6\\  are  consider- 
ably more  jagged  than  the  configuration  space  perturbations,  a  phenomenon 
long  observed  for  the  evolution  of  the  total  3iV-dimensional  perturbations  Ar 
and  Av  [26].  However,  a  comparison  of  the  individual  |<$r|'s  and  |£v|'s  yields 
qualitative  differences  for  different  values  of  e.  For  the  smallest  softening  pa- 
rameter, namely  e  =  0.005,  the  kicks  in  the  velocity  perturbation  tend  to  occur 
at  the  same  times  as  the  kicks  in  the  configuration  space  perturbation.  How- 
ever, for  e  =  0.18,  the  kicks  in  velocity  tend  to  be  completely  out  of  phase  with 
the  kicks  in  position.  The  delay  in  the  velocity  perturbation  of  the  kicks  due 
to  short  range  encounters  becomes  significant  when  there  is  a  large  amount  of 
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softening.  For  e  =  0.18,  it  was  also  observed  that  the  "bumps"  in  the  configura- 
tion and  velocity  space  perturbations  sometimes  appear  oscillatory:  some  15% 
of  the  randomly  chosen  particles  can  be  said  to  display  this  sort  of  behavior. 

Since,  for  such  a  large  softening,  the  perturbation  never  grows  very  large 
(the  typical  \6r\  <  10-6  even  at  the  end  of  the  simulations),  it  seemed  possible 
that  these  "oscillations"  might  be  associated  with  perturbations  of  some  binary 
system.  This  hypothesis  was  in  fact  studied  in  some  detail,  but  the  results  were 
inconclusive.  It  was  found  that,  although  a  single  other  particle  did  often  tend 
to  remain  in  the  neighborhood  of  the  "oscillating"  particle  during  the  entire 
run,  the  times  at  which  this  particle  approached  the  first  could  not  be  simply 
correlated  with  the  kicks  in  the  perturbation  of  the  original  particle.  This 
problem  clearly  warrants  future  attention. 

Overall,  the  typical  growth  of  individual  perturbations  is  at  least  roughly 
fit  by  a  simple  exponential  form  |6r(tf)|  ~  |6r(0)|  exp(t/t\).  However,  as  for 
the  total  3iV-dimensional  perturbation,  the  growth  time  t\  clearly  depends 
on  the  value  of  c.  When  probed  over  a  time  of  approximately  2.8tcr,  the 
mean  growth  times  t\  are  comparable  to,  but  slightly  longer  than,  the  growth 
times  t\  for  the  full  3iV-dimensional  perturbations.  Indeed,  one  finds  that, 
for  e  =  0.005,  tx  =  0.17±0.01*cr,  for  e  =  0.05  that  tx  =  0.59*cr,  and,  for 
e  =  0.18  that  t\  =  1.41±0.37£cr.  These  values  are  to  be  compared  with  the 
time  scales  tA>L  =  0.15±0.01*cr,  0.52±0.02<cr,  and  1.17±0.21<cr  and  the  corre- 
sponding t\jj  recorded  in  Table  2.  The  fact  that  these  individual  time  scales 
are  slightly  longer  reflects  the  fact  that  the  full  3iV-dimensional  perturbation  is 
dominated  by  those  particles  whose  perturbation  \8r\  is  growing  more  quickly 
than  average.  When  probed  only  over  the  initial  0.7icr,  the  corresponding  time 
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scales  t\  are  increased  by  roughly  a  factor  of  two.  In  this  case,  one  finds  respec- 
tively the  values  tx  =  0.23±0.04<cr,  1.01<cr,  and  2.40±0.09*cr.  This  difference 
reflects  the  fact  that,  initially,  the  perturbation  either  decreases  or,  at  best, 
increases  only  very  slowly. 

A  related  question  concerns  the  degree  to  which  the  perturbations  of  dif- 
ferent particles  in  fact  grow  at  roughly  the  same  rate.  For  e  =  0.005,  the 
growth  rates  tend  to  agree  to  within  about  10%  on  a  time  scale  of  2.8tfCr  and 
to  within  ~  50%  over  0.7tcr,  but,  over  the  first  0.35rcr,  only  at  the  175%  level. 
For  e  =  0.18,  the  agreement  is  at  the  10%  level  over  5.6£cr,  at  the  50%  level 
over  2.8*cr,  and  at  the  350%  level  over  0.7tcr.  Over  the  initial  0.35tfCr,  the 
agreement  actually  improves  to  ~  75%. 

This  agreement  during  the  initial  transient  period  for  e  =  0.18  is  related 
to  the  fact  that,  for  a  majority  of  the  particles,  the  perturbation  actually 
decreases  initially.  This  decrease  is  in  fact  observed  for  every  value  of  e,  but 
it  becomes  more  significant  for  the  larger  values:  the  initial  decrease  is  almost 
always  larger  for  e  =  0.18,  and  the  frequency  of  its  occurrence  is  greater. 
Moreover,  the  perturbation  reaches  its  minimum  earlier  for  the  lowest  values 
of  e.  For  e  =  0.005,  the  minimum  almost  always  occurs  before  0.35rcr  whereas, 
for  e  =  0.18,  this  usually  occurs  between  0.35  and  0.7tcr.  For  e  =  0.005,  the 
initial  growth  is  negative  for  some  60%  of  the  particles,  for  the  0.05  case,  ~  70% 
of  the  particles,  and,  for  e  =  0.18,  ~  80%  of  the  particles. 

For  those  particles  for  which  the  perturbation  does  not  initially  decrease, 
there  is  almost  invariably  only  a  smooth,  very  slow  initial  growth,  correspond- 
ing typically  to  a  time  scale  ~  3  -  btCT.  There  are  almost  no  instances  in  which 
the  perturbation  begins  to  grow  very  quickly  at  the  outset.   Rapid  growth  in 
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the  perturbation  appears  to  begin  rather  abruptly  at  some  later  time,  pre- 
sumably after  an  initial  kick  associated  with  some  close  encounter.  This  fact 
reinforces  the  picture  that  both  close  encounters  and  the  bulk  mean  field  play 
an  important  role  in  the  instability.  It  is  evident  that,  without  close  encoun- 
ters, the  perturbation  will  only  grow  very  slowly,  but  it  is  also  clear  that  the 
mean  field  plays  a  role  in  amplifying  the  effects  of  the  close  encounters. 


CHAPTER  4 

COLLISIONAL  RELAXATION  IN  A 

NONINTEGRABLE  POTENTIAL 

It  was  seen  in  the  previous  chapter  that  the  instability  of  the  JV-body 
gravitational  problem  depends  on  both  the  mean  field  and  'noise'  due  to  close 
encounters,  independent  of  the  size  of  the  system,  close  encounters  remain  im- 
portant even  when  N  is  large,  and  decrease  the  time  scale  for  instability.  This 
would  imply  that  for  galactic  systems,  the  collisionless  Boltzmann  equation, 
as  discussed  in  section  1.2,  may  not  be  applicable  when  determining  the  time 
scale  for  instability. 

Alternatively,  one  could  formulate  a  Langevin  equation,  which  incorporates 
the  effects  of  close  encounter  noise  directly.  The  standard  pardadigm  is  to 
model  these  close  encounters  essentially  as  a  Brownian  process,  while  ignoring 
the  effects  of  the  bulk  mean  field  [27,28].  This  decoupling  of  the  mean  field  and 
the  fluctuations  due  to  close  encounters  implies  that  the  collisional  relaxation, 
which  was  discussed  in  section  1.1,  is  much  greater  than  the  age  of  the  universe 
for  galactic  systems. 

This  conclusion  is  a  direct  consequence  of  the  assumption  that  fluctuations 
due  to  close  encounters  can  be  completely  decoupled  from  the  bulk  mean  field. 
This  may  indeed  be  true  for  systems  in  which  the  trajectories  are  stable.  Non- 
integrable  systems,  however,  are  characterized  by  the  existence  of  exponentially 
unstable  trajectories.  Thus,  the  standard  calculation  of  the  relaxation  time  of 
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a  particle  in  a  mean  field  is  reexamined  here,  assuming  the  mean  field  tra- 
jectories are  exponentially  unstable,  and  characterized  by  a  positive  Liapunov 
exponent. 

4.1  Relaxation  in  an  Inverted  Harmonic  Oscillator 
Here  we  consider  collisional  relaxation  in  an  exactly  soluble  toy  model, 
the  inverted  harmonic  oscillator,  in  which  all  trajectories  are  exponentially 
unstable.   The  calculation  is  similar  to  that  of  Chandrasekhar    [27],  and  can 
be  extracted  from  this  via  an  analytic  continuation  u2  — ♦  —  u>  . 
Consider  the  generalized  one- dimensional  Langevin  equation 

d  x         dx         o  ,  v  , 

_  +  „__w2  +  F(()  =  0,  (4.1) 

where  F(t)  is  the  stochastic  force  and  r\  is  the  associated  dynamical  friction.  As 
a  simplifying  idealization,  assume  that  F(t)  is  delta-correlated  white  noise  with 
zero  mean,  so  that  F  and  rj  are  related  via  a  fluctuation-dissipation  theorem 

(F(<1)F(i3))=e^(*i-<2),  (4.2) 

where  0  defines  a  characteristic  'temperature',  i.e.  a  typical  mean  squared  ve- 
locity. Also  assume  that  r),  the  coefficient  of  dynamical  friction,  is  independent 
of  position  and  velocity. 

It  is  straightforward  to  right  down  the  solutions  to  an  ordinary  differen- 
tial equation  for  the  damped  inverted  oscillator,  with  F(t)  being  considered  a 
driving  term,  ignoring  for  the  moment  the  fact  that  F(t)  is  a  stochastic  force. 
The  solutions  to  the  homogeneous  equation  can  be  expressed  as 

x  =  aie"1'  +  a2e"2<  (4.3) 
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where  /xi  and  ^2  are  the  roots  of  the  characteristic  equation 

H2  +  rjn  -  u2  =  0.  (4.4) 

To  determine  a  particular  solution,  it  is  convenient  to  use  the  method  of 
variation  of  parameters,  as  discussed  in,  for  example,  [29].  Thus  a\  and  02  are 
replaced  by  functions  of  time  and  determined  by  the  equations 

(He"1*  +  a2e"2t  =  0  (4.5) 

and 

/Jidie"1'  +  /i2a2e/i2<  =  F(t).  (4.6) 

The  equations  (4.5)  and  (4.6)  can  be  used  to  obtain  integral  solutions  for 
a\(t)  and  02(2),  which  can  be  inserted  into  equation  (4.3)  to  obtain 

x  -  x0e-«</2  cosh  £  -  £°2  +  2!V*/2  sinh  &  m  f  d(F{0m         (4.7) 
2  P  2       J0 

and 

v  -  voe-**'2  cosh  £  +  XZ*!^-*fl  sinh  §  .  /'  deF(0*(0,      (4-8) 

where  /32  =  rj2  +  4u>2.  The  right-hand  sides  of  equations  (4.7)  and  (4.8)  involve 
the  functions 

*(0  =  y-^-W  sinh[/9(t  -  0/2]  (4.9) 

and 

$  =  rf¥/d*.  (4.10) 

At  this  point,  is  important  to  consider  that  F(t)  is  a  stochastic  force,  and 
hence  the  equations  (4.7)  and  (4.8)  can  only  be  solved  in  the  sense  of  specifying 
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a  probability  distribution  function  which  can  be  used  to  extract  average  quan- 
tities for  the  variables  of  interest.  Given  the  form  of  the  left-hand  sides  of  equa- 
tions (4.7)  and  (4.8),  the  distributions  functions  W(x,t;xo,vo),W(v,t;xo,vo), 
and  W(x,  u,  t;  xq,vq)  can  be  written  in  terms  of  the  integrals 

/  d£*2(0;     /  de*2(0   and     /  df*(0*(0  (4.11) 

Jo  Jo  Jo 

as  proven,  for  example,  by  Chandrasekhar    [27]. 

More  straigtforwardly,  the  averages  of  the  moments  (xpv9)  can  be  computed 
directly  from  equations  (4.7)  and  (4.8).  Since  the  stochastic  force  has  zero 
mean,  averaging  the  equations  (4.7)  and  (4.8)  immediately  gives 

x  =  zee""'/*  cosh  f  -  m  +  **e-«/*  sinh  |  (4.12) 

and 

v  =  voe-*l>  cosh  &  +  ^-2"2*0e-*/2  sinh  #  (4  13) 

while  the  second  moments  can  be  obtained  using  the  integrals  (4.11). 

In  an  astronomical  context,  the  instability  would  be  expected  to  grow  on 
the  time  scale  of  a  fiducial  crossing  time,  i.e.  the  dynamical  time  tj),  so  u 
would  be  roughly  comparable  to  tp.  The  coefficient  of  dynamical  friction, 
which  is  related  to  the  noise  via  a  fluctuation-dissipation  theorem,  would  be 
determined  by  the  time  scale  set  by  close  encounters,  so  that,  by  definition, 
•q  =  t~^  ,  where  tji  is  the  classical  binary  relaxation  time.  Recall  from  chapter 
1  that  tR  is  much  greater  than  the  age  of  the  universe,  so  that  tj)  <Ctf.fi. 

It  is  also  useful  to  neglect  transients  and  focus  on  the  behavior  that  occurs 
after  several  crossing  times,  as  well  as  neglecting  effects  which  proceed  on  the 
much  longer  time  tR.  Thus,  we  wish  to  assume  that  u>-1  <  t  <  r/_1. 
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In  this  approximation,  the  noise  and  friction  are  weak  and  fi  «  2uj.  Since 
we  are  neglecting  transients,  sinhwt  ss  ewt/2.   Thus,  the  solutions  (4.12)  and 
(4.13)  reduce  to 

(x)  =  \  exp[(u,  -  r,/2)t]  (*„  +  £  +  ^x0)  (4.14) 


and 


(v)  =  -  exp[(u;  -  rj/2)t]  (vq  +  ux0  -  ^-"oj  •  (4.15) 


This  result  reduces  to  the  deterministic  result  F(t)  — ►  0  when  the  friction 
coefficient  rj  <C  2u,  so  that  2u  —  tj  w  2u>. 

In  the  weak  noise  and  weak  friction  limit  /?  w  2u>,  the  integrals  (4.11) 
simphfy  as  well  and  can  be  written  as 

J  d&2(0  a  ^  exp[(2u;  -  //)<],  (4.16) 


/ 


H/-        1 


^*(0«£jexp[(2u; ->;)<],  (4-17) 
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and 

jf  de*(0*(0  «  8^2  exp[(2W  -  f,)i).  (4.18) 

With  the  help  of  the  fluctuation-dissipation  theorem  (equation  (4.2)),  we  can 
then  write  the  second  moments  as 

<*2>  =  <*>2  +  £^exp[(2o,  -  r/)<],  (4.19) 

(„2)  =  {vf  +  ^exp[(2u;  -  r,)t],  (4.20) 
and 

(xv)  =  (x)  (v)  +  ^-Xexp[(2a;  -  Vt)].  (4.21) 

These  moment  equations  can  now  be  used  to  determine  approximations  for 
the  growth  rates  for  the  mean  squared  changes  in  position  and  velocity  (\8x\2) 
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and  (|£t>|2),  as  well  as  (SxSv).   To  lowest  order  in  rj,  {\Sx\2)  «  (x2)  —  (x)2  so 
that,  for  2u>  —  rj  «  2a;, 

<l**|2)  =  |J£e™  (4.22) 

and  similarly 

(l*vla)  =  §V*  (4.23) 

and 

<*s*y>  =  Ip"".  (4.24) 

Recall  from  chapter  1  that  the  relaxation  time  was  defined  in  terms  of 
changes  in  the  mean  square  velocity,  the  classical  result  being 

JftCl  ~  ^LL.  (4.25) 

This  can  be  contrasted  to  the  current  estimate  obtained  from  equation  (4.23), 
using  the  facts  that  the  'temperature'  0  corresponds  to  a  mean  square  velocity 
(v2),  u  ~  t~p,  rf  =  <£*,  and,  from  the  virial  theorem,  |*  ~  g^  (equation  1.6). 
The  conclusion  is  that  collisional  effects  in  an  inverted  harmonic  oscillator,  for 
which  the  mean  field  orbits  are  unstable,  induce  changes  in  the  mean  square 
velocity  of  trajectories  which  are  estimated  to  be 

<H!)^exp(*/(o).  (4.26) 

Similar  conclusions  are  obtained  for  changes  in  the  other  moments  as  well. 

It  is  apparent  that,  even  if  N  is  very  large,  the  exponential  factor  in  equa- 
tion (4.26)  implies  that  the  time  scale  for  appreciable  changes  in  velocity  will 
be  relatively  short.  Thus,  the  effective  relaxation  time  for  trajectories  in  a 
nonintegrable  potential  can  be  orders  of  magnitude  smaller  than  the  classical 
binary  relaxation  time,  for  which  stable  mean  field  orbits  are  implied. 


55 
The  above  conclusion  is  not  necessary  valid  for  changes  in  energy,  however. 
The  mean  change  in  particle  energy  can  be  estimated  from  the  above  moment 
equations  to  yeild,  for  (6E)  =  \{v2  -u2x2)  -  ^(u^  -  w2x^y), 

(M)  =  -{^)(vl-»2*iywi-  (4-27) 

For  this  toy  model,  perturbation  in  particle  energy  induced  by  the  noise  and 
friction  will  indeed  grow  exponentially.  However,  this  is  a  consequence  of  the 
fact  that,  for  this  model,  the  unperturbed  positions  and  velocities  diverge  for 
late  times,  even  though  the  unperturbed  particle  energy  is  in  fact  conserved,  for 
a  generic  potential,  this  is  not  true  and  the  natural  time-scale  associated  with 
changes  in  particle  energy  is,  in  fact,  the  classical  relaxation  time,  tji  =  17     . 

While  the  above  conclusions  were  arrived  at  for  the  case  of  an  exactly 
soluble  toy  model,  the  same  basic  idea  can  be  exploited  for  any  generic  potential 
with  unstable  mean  field  orbits.  This  is  the  idea  explored  in  the  next  section. 

4.2  Relaxation  in  a  Generic  Nonintegrable  Potential 
Consider  a  particle  moving  in  a  potential  U(r)  so  that  it  follows  a  smooth, 
deterministic  trajectory  ra(t)  which  obeys 

d2ra       dU(r) 

-w  +  -^-  =  0   (fl  =  1'2'3)-  (4-28) 

The  background  potential  U(r)  may  be  nonintegrable  and  have  some  exponen- 
tially unstable  orbits.  If  we  now  include  a  fluctuating  force  and  its  associated 
friction,  we  obtain  a  generalized  Langevin-type  equation  of  the  form 

d2Ra       dU(R)  dRa       „,  ,        , 

solutions  of  which  represent  noisy  trajectories  Ra(t). 
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The  primary  interest  is  in  the  difference  between  these  trajectories,  8ra  = 
Ra  —  r°.  Treating  6ra  as  a  perturbation  and  linearizing,  we  obtain 

#6r*         6r*      #U(r)     b 

-W  +  V-dT  +  -dr^8T   ~F   'VV  '  (4-30) 

where  va  =  dra/dt. 

Since  we  are  assuming  the  orbit  ra  is  exponentially  unstable,  at  least  one 
of  the  eigenvalues  of  the  matrix  ~  t £1  must  be  negative,  and  hence  of  the  form 
—u!2  where  u>2  >  0.  It  is  reasonable  to  assume  that  the  stochastic  force  Fa 
acts  randomly  in  all  directions,  so  that  it  will  eventually  probe  the  direction 
associated  with  u2A.  In  any  event,  if  the  orbit  is  exponentially  unstable  in  all 
directions,  so  that  for  each  value  of  a  the  eigenvalues  are  of  the  form  —  u>2,  con- 
sidering only  the  smallest  value  of  u>2  will  only  underestimate  the  importance 
of  the  instability.  Thus,  it  is  reasonable  to  replace  »  h^T'a  by  —  oj28a  ,  where  u> 
can  be  considered  the  smallest  value  of  u\. 

The  value  u  can  be  interpreted  as  a  Liapunov  exponent  for  the  mean  field 
trajectory,  which  must  be  positive  for  exponentially  unstable  directions.  While 
we  are  interpreting  uj  to  be  the  smallest  Liapunov  exponent  to  be  safe,  it 
should  be  apparent  that  the  direction  associated  with  the  largest  Liapunov 
exponent  should  quickly  dominate  the  instability,  and  hence  should  be  the 
physically  relevant  direction.  Indeed,  numerical  simulations  indicate  that  the 
largest  Liapunov  exponent  generally  provides  a  reasonable  characterization  of 
the  instability,  as  shown  in,  e.g.     [30]. 

The  linearized  equation  can  now  be  written 

_r+,__u^=F._v,..  (4.S1) 
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which  is  an  inhomogeneous  stochastic  differential  equation  for  Sra,  with  the 
quantity  Fa—r)va  acting  as  the  source.  This  equation  is  very  similar  to  equation 
(4.1)  for  the  toy  model.  Thus  it  is  straightforward  to  write  down  solutions  to 
equation  (4.31)  in  terms  of  the  functions  ^  and  $  given  in  equations  (4.9)  and 
(4.10).  The  resulting  equations  are 

sr*=  /'deiF-co-^'cowo  (4-32) 

Jo 


and 

d6r° 


dt    =8v°  =  j    d£[Fa(0-Vva(0m)-  (4-33) 

In  order  to  evaluate  the  moments  of  this  equation,  it  is  again  convenient 
to  assume  delta-correlated  Gaussian  white  noise  with  zero  mean,  so  that  the 
three-dimensional  fluctuation-dissipation  theorem  becomes 

(F(h)-F(t2))  =  erjSD(t2-ti).  (4.34) 

Therefore  equation  (4.32)  can  be  used  to  determine  the  second  moment 

<|*r|2>  .  (  j  di  J   dC*(0*(C)[**(0  ~  Wa(0)[Fa(0  ~  VVa(0})       (4-35) 

Since  the  friction  and  the  noise  are  uncorrelated  and  the  noise  is  assumed  to 
have  zero  mean,  this  can  be  written 

(\sr\2)  =  ev  I  d&2(o  +  A  I  dev(o*(0  2-  (4-36) 

Jo  i  Jo 

Similar  results  can  of  course  be  obtained  for  (|<5v|2)  and  {6r  •  <5v). 

It  is  appropriate  to  evaluate  these  quantities  in  the  late  time,  weak  noise 
limit  considered  in  the  last  section.  In  this  case,  the  second  term  in  equation 
(4.36)  can  be  neglected  for  typical  values  of  the  velocity,  since  it  is  of  order  t/2. 
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Since  this  term  is  positive  definite,  by  doing  so  we  axe  only  underestimating 
the  growth  rate.  The  remaining  integral  is  given  by  equation  (4.16)  and  hence 
the  result  is 


Similarly, 


<l*r|2)  =  |^.  (4.37) 


<l*v|2)  =  fV"".  (4-37) 


and 

(6r.6v)  =  ^e™.  (4.37) 

We  thus  recover  the  same  results  as  for  our  model  in  section  2,  namely,  that 
the  relaxation  time  associated  with  changes  in  the  mean  square  of  the  phase 
space  variables  may  be  orders  of  magnitude  shorter  than  the  classical  result. 


CHAPTER  5 
NONHYPERBOLIC  SHADOWING  IN  MAPPINGS 

The  fact  that  chaotic  dynamical  systems  exhibit  sensitive  dependence  on 
initial  conditions  implies  that  the  presence  of  even  weak  noise  may  radically  al- 
ter the  details  of  individual  orbits.  We  have  seen  in  previous  chapters  (chapter 
3  and  chapter  4)  how  noise  effects  the  time  scale  of  the  instability  in  a  nonin- 
tegrable  potential.  But  how  does  the  noise  affect  the  details  of  the  individual 
orbits?  More  importantly,  even  though  the  individual  orbits  are  drastically 
altered  by  the  presence  of  noise,  are  the  macroscopic  properties  of  the  system 
altered  as  well?  In  other  words,  do  the  noisy  orbits  behave  the  same  as  the 
deterministic  orbits,  in  a  statistical  sense,  or  is  the  behavior  of  the  orbits  in  a 
noisy  chaotic  system  qualitatively  changed? 

The  answer  to  these  questions  is  of  extreme  importance:  if  the  behavior  of 
the  noisy  orbits  is  qualitatively  similar  to  the  orbits  in  the  deterministic  system, 
we  might  expect  that  the  statistics  of  the  orbits  obtained  for  both  systems 
should  be  identical.  And  since  we  know  that  for  a  chaotic  system  the  character 
of  the  individual  orbits  are  essentially  useless,  due  to  the  exponential  divergence 
of  trajectories,  the  addition  of  noise  would  not  hinder  any  understanding  of  the 
orbits  in  this  regard.  For  such  systems,  why  bother  analysing  the  noisy  system 
at  all,  since  one  would  expect  to  obtain  identical  results  for  the  corresponding 
deterministic  system. 
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The  above  scenario  is  not  always  true  however,  and  it  is  important  to 
understand  the  qualitative  change  that  noise  has  on  typical  orbits  in  a  non- 
integrable  system.    In  the  next  few  chapters,  I  will  attempt  to  answer  these 
questions  using  shadowing  theorems. 

5.1  Shadowing  as  Noise  Reduction 
In  many  cases,  noise  is  an  unwanted  and  annoying  artifact  in  physics,  and 
it  is  often  hoped  that  one  can  minimize  its  effects.  Undesirable  noise  may 
be  introduced  into  physical  experiments  via  errors  associated  with  the  mea- 
surement process,  and  into  computer  experiments  via  roundoff  and  truncation 
errors.  Noise  due  to  measurement  errors  can  generally  be  considered  obser- 
vational noise,  sometimes  called  additive  noise,  in  which  noise  introduced  at 
one  instant  does  not  affect  the  orbit  at  the  next  instant.  In  this  case,  given  a 
dynamical  system  xn+i  =  /(x„),  a  noisy  orbit  is  related  to  the  deterministic 
orbit  by  simply  adding  a  noise  term  tj,  pn  =  xn  +  rjn. 

f  f1 


t-1 

Figure  5.1  Noisy  measurements  made  at  three  different  times.  When  these 
measurements  are  transported  to  the  same  point  in  time,  they  distort  in  ac- 
cordance with  the  local  derivatives  of  the  dynamical  system. 
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In  this  case,  there  is  a  straightforward  method  to  reduce  unwanted  noise 
called  nonlinear  smoothing.  Consider  making  noisy  measurements  at  times 
t-1,  t,  and  t+1,  as  described  by  Farmer  and  Sidorowich  [31],  and  shown  in 
figure  5.1.  The  measurement  at  time  t-1  is  evolved  forward  in  time  to  time 
t,  and  expands  along  the  unstable  manifold  while  contracting  along  the  stable 
manifold.  Evolving  the  measurement  at  time  t+1  backwards  to  time  t  causes 
the  data  to  expand  along  the  stable  manifold  and  contract  along  the  unstable 
manifold.  The  intersection  of  these  three  data  sets  at  time  t  yield  a  more  refined 
measurement.  The  amount  that  the  transported  data  set  distorts  depends  on 
the  rate  of  contraction  or  expansion  along  the  stable  or  unstable  manifold, 
taking  into  account  the  fact  that,  for  Hamiltonian  systems,  the  area  of  the 
data  region  must  remain  the  same.  This  scheme  fails  when  these  manifolds  do 
not  lie  in  distinct  directions,  a  difficulty  which  will  be  discussed  more  later. 

By  contrast,  Dynamical  noise  introduced  at  one  instant  alters  the  dynamics 
of  the  orbit  in  the  future,  and  in  a  chaotic  system  causes  a  divergence  from  the 
purely  deterministic  orbit. Computer  errors  are  often  the  cause  of  dynamical 
noise.  Roundoff  errors  are  due  to  the  finite  floating  point  arithmetic  which 
occurs  whenever  one  does  a  computer  simulation  involving  irrational  numbers, 
whereas  truncation  errors  occur  whenever  attempting  to  integrate  a  differential 
equation  numerically,  and  depend  on  the  details  of  the  integration  scheme  used. 

Even  so,  the  effect  of  computer  errors  on  the  behavior  of  particle  orbits 
is  often  ignored  by  the  galactic  dynamicist  doing  numerical  simulations.  It  is 
frequently  assumed  that,  even  if  the  details  of  the  noisy  orbits  differ  drastically 
from  the  deterministic  orbits,  the  properties  of  the  orbits  are  the  same,  in  some 
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statistical  sense.  Statistical  tests  seem  to  indicate  that  this  is  indeed  true  for 
computer  simulations  of  N-body  gravitational  systems    [32,33]. 

Moreover,  in  some  cases  the  neglect  of  computer  errors  in  a  simulation  can 
perhaps  be  justified  by  the  possibility  that  they  may  actually  mimic  the  random 
perturbations  that  influence  real,  physical  systems.  In  the  case  of  galactic  dy- 
namics, these  perturbations  may  be  assumed  to  be  due  to  interactions  amongst 
the  various  components  of  the  galaxy,  for  example,  close  encounters  of  nearby 
stars.  While  it  may  be  true  that  the  computer  noise  will  not  exactly  reproduce 
the  physical  noise,  both  types  of  noise  are  dynamical,  and  the  exact  form  of 
the  noise  may  not  really  matter.  If  this  be  true,  one  should  not  necessarily  be 
concerned  with  the  use  of  noise  reduction  methods  to  remove  the  impact  of 
noise  on  a  simulation.  Instead,  the  concern  should  be  whether  the  amplitude 
of  the  noise  has  any  physical  relevance.  The  existence  of  computer  errors,  or 
any  dynamical  noise  for  that  matter,  should  not  necessarily  be  considered  de- 
structive in  a  computer  simulation,  but  may  actually  be  useful  for  modelling  a 
given  physical  system. 

Nevertheless,  the  galactic  dynamicist  should  still  be  concerned  as  to  the 
qualitative  effect  a  given  amount  of  noise  will  have  on  the  properties  of  orbits 
in  a  chaotic  dynamical  system.  A  potentially  useful  technique  for  determining 
the  qualitative  effect  of  noise  on  a  dynamical  system  may  be  to  study  the 
shadowability  of  the  system.  Shadowable  systems  possess  the  property  that, 
in  a  certain  sense,  a  deterministic  orbit  remains  close  to  a  noisy  orbit.  What 
exactly  is  meant  by  the  term  shadowability  will  be  made  clear  in  the  next 
sections. 


63 
While  shadowing  was  originally  motivated  by  the  need  for  a  noise  reduc- 
tion scheme  for  chaotic  systems,  I  will  be  interpreting  the  shadowing  property 
instead  as  a  measure  of  the  change  in  qualitative  behavior  of  trajectories  in  the 
noisy  system,  as  opposed  to  the  deterministic  system.  In  other  words,  whereas 
the  traditional  discussion  of  shadowing  considered  the  deterministic  orbit  to 
be  the  true  orbit  of  the  system  of  interest  while  the  noisy  orbit  is  considered 
the  undesired  orbit,  I  will  be  considering  the  noisy  orbit  to  be  the  orbit  of  in- 
terest, with  the  goal  of  determining  how  well  a  low-dimensional,  deterministic 
system  can  model  the  more  relevant,  noisy  physical  system.  In  what  follows,  I 
will  be  using  the  traditional  definitions  for  the  orbits,  so  readers  familiar  with 
shadowing  arguments  of  other  authors  will  be  able  to  follow  the  discussion. 

5.2  Hyperbolic  Shadowing 
Hyperbolic  systems,  as  mentioned  in  Chapter  Two,  are  systems  for  which 
the  stable  and  unstable  manifolds  always  intersect  transversely.  Such  systems 
have  been  proven  to  be  shadowable  for  all  time  by  Anosov    [34]  and  Bowen 
[35]. 

Theorem.  (Anosov-Bowen,  1967)  For  hyperbolic  systems,  given  a  shadow 
distance  e  >  0,  there  exists  a  8  >  0  such  that  any  6 -pseudo-orbit  can  be 
e-shadowed  by  a  true  orbit  for  all  time. 

This  theorem  concerns  the  orbits  of  a  hyperbolic  map  xn+i  =  /(xn).  A 
true  orbit  is  defined  to  be  the  set  of  points  {xn}£L0,  such  that  x„+i  =  /(xn) 
for  n  =  0, . . . ,  JV  —  1.  Likewise,  a  6-pseudo-orbit  is  the  set  of  points  {yn}n=o 
such  that  |yn+i  -  f(yn)|  <  6  for  n  =  0, . . .  ,JV  -  1.  In  other  words,  the  6- 
pseudo-orbit  is  the  noisy  orbit  with  noise  amplitude  6;  it  remains  within  6  of 
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shadow  orbit 


deterministic  orbit 
Figure     5.2     Deterministic  and  noisy  orbits  of  a  dynamical  system.    The 
lowest  curve  is  the  purely  deterministic  orbit  for  a  given  initial  condition.  The 
middle  curve  is  the  noisy  orbit  with  the  same  initial  condition.    The  upper 
curve  is  the  shadow  orbit  that  remains  within  a  distance  e  of  the  noisy  orbit. 


the  true  orbit  which  the  system  would  have  followed  at  that  time  step  without 
any  noise. 

Now,  the  e-shadow  orbit  is  a  true  orbit  of  x„+i  =  /(xn)  if  |yn  —  xn|  <  e 
for  n  =  0,  ...,JV  —  1.  The  distance  e,  the  maximum  distance  between  the 
shadow  orbit  and  the  pseudo-orbit,  will  be  called  the  shadow  distance.  The 
hyperbolic  shadowing  theorem  states  that  a  deterministic  orbit  and  a  noisy 
orbit  remain  within  a  distance  e  as  long  as  they  are  separated  by  no  more  than 
6  initially.  Thus,  the  shadow  orbit  and  the  noisy  orbit  do  not,  in  general,  have 
the  same  initial  conditions,  and,  as  such,  the  initial  conditions  of  an  orbit  in 
a  hyperbolic  dynamical  system  need  not  be  known  exactly  to  determine  the 
qualitative  behavior  of  the  orbit.  This  is  shown  in  figure  5.2. 

It  has  been  shown  by  Anosov  [36]  that  hyperbolic  systems  are  structurally 
stable,  i.e.  a  perturbed  hyperbolic  system  remains  hyperbolic.  For  such  sys- 
tems, one  would  not  expect  any  change  in  the  statistical  quantities  calculated 
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for  the  original  hyperbolic  system  and  the  perturbed  system.  It  is  thus  ap- 
parent that  the  addition  of  noise  in  a  hyperbolic  system  will  not  induce  any 
qualitative  change  in  the  behavior  of  orbits  in  the  system,  even  though  the 
noisy  orbits  are  drastically  different  from  the  deterministic  ones.  This  is  true 
irrespective  of  the  form  of  the  noise:  observational  or  dynamical,  white  or 
colored,  gaussian  or  not. 

The  shadowing  theorem  obviously  holds  for  any  system  with  observational 
noise,  hyperbolic  or  not.  In  this  case,  since  the  presence  of  the  noise  does 
not  alter  the  dynamics,  the  shadow  distance  is  then  merely  equal  to  the  noise 
amplitude.  The  power  in  the  shadowing  theorem  lies  in  its  predictions  for  the 
case  of  dynamical  noise. 

It  is  important  to  note  that,  in  the  standard  interpretation,  the  true  orbit 
in  the  Anosov-Bowen  theorem  is  only  a  true  orbit  of  the  purely  deterministic 
dynamical  system.  However,  if  one  wishes  to  consider  the  noise  to  be  of  some 
physical  origin,  then  the  ^-pseudo-orbit  should  more  appropriately  be  consid- 
ered a  true  orbit  of  the  noisy  dynamical  system.  Given  this  new  interpretation, 
the  hyperbolic  shadowing  theorem  states  that  a  true,  noisy  orbit  with  noise 
amplitude  8  would  remain  within  e  of  some  orbit  of  the  purely  deterministic 
system. 

5.3  Shadowing  of  Nonhvperbolic  Systems 

However,  most  problems  of  physical  interest  are  not  hyperbolic,  and  the 

Anosov-Bowen  theorem  is  not  applicable  to  the  majority  of  systems  that  the 

physicist  or  galactic  dynamicist  would  like  to  represent.   Two  representative 

systems  which  are  considered  paradigms  of  chaotic  systems  are  the  Henon 
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map  and  the  standard  map,  neither  of  which  is  hyperbolic.  Other  com- 
monly studied  systems  of  physical  interest  which  are  nonhyperbolic  are  the 
Henon-Heiles  system,  often  used  as  a  simple  galactic  model,  and  the  complex 
Ginsburg- Landau  equation,  which  has  important  applications  in  condensed 
matter  physics.  Hence,  a  more  general  shadowing  theorem  is  necessary  to 
determine  the  shadowability  of  nonhyperbolic  systems. 

The  difficulty  in  shadowing  nonhyperbolic  systems  lies  in  the  existence 
of  homoclinic  tangencies  and  near-tangencies.  Homoclinic  tangencies  occur 
when  the  angle  between  the  stable  and  unstable  manifolds  of  a  point  along  a 
trajectory  approaches  zero.  Referring  back  to  figure  1,  this  means  that  the 
ovals  representing  the  noisy  measurements  lie  along  the  same  axis.  How  close 
to  a  homoclinic  tangency  a  trajectory  is  able  to  get  depends  strongly  on  the 
noise  amplitude.  Referring  to  figure  3,  as  the  angle  6„  approaches  zero,  the 
trajectory  approaches  a  homoclinic  tangency.  If  9n  is  small  enough  so  that  the 
noise  allows  the  trajectory  to  "jump"  over  a  stable  manifold,  it  will  then  diverge 
from  the  deterministic  trajectory.  The  system  has  approached  a  near-tangency, 
and  the  trajectory  is  no  longer  shadowable. 

Thus,  the  trajectories  in  nonhyperbolic  systems  cannot  be  shadowed  for 
all  time,  and  there  is  no  a  priori  gaurantee  that  the  noisy  and  deterministic 
trajectories  behave  in  a  similar  fashion.  However,  trajectories  can  be  shadowed 
as  long  as  they  remain  in  regions  which  are  "hyperbolic  enough,"  that  is,  as 
long  as  they  remain  far  enough  away  from  homoclinic  near-tangencies.  The 
length  of  time  an  orbit  remains  in  such  a  sufficiently  hyperbolic  region  will  be 
called  the  shadowing  time.  The  'shadowing  time'  may  be  defined  as  the  length 
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of  time  N  such  that  a  shadow  orbit  exists  given  a  noisy  orbit  for  all  times 
n  =  0, 1, ...,  N,  but  such  an  orbit  does  not  exist  for  n  =  0, 1, ...,  iV  +  1. 

Several  algorithms  have  been  implemented  that  calculate  shadowing  times 
for  one-  and  two-dimensional  maps.  In  general,  these  procedures  rely  on  a 
refinement  procedure  which  will  be  described  in  the  next  section.  These  algo- 
rithms have  limitations  that  are  avoided  in  a  more  powerful  algorithm  due  to 
Sauer  and  Yorke  [37].  They  derive  a  nonhyperbolic  shadowing  theorem  that 
suggests  an  algorithm  which  can  be  used  to  calculate  how  long  an  orbit  can  be 
shadowed  before  approaching  a  near-tangency,  for  a  general  map  /  of  arbitrary 
dimension. 

The  Sauer  and  Yorke  algorithm  relies  on  determining  how  small  the  angle 
between  the  stable  and  unstable  manifolds,  6n,  is  allowed  to  get  before  a  near- 
tangency  is  encountered.  This  requires  information  regarding  the  linearized 
map;  specifically,  the  bounds  on  the  growth  of  /„  in  the  direction  S„  and  of  /"* 
in  the  direction  Un.  For  this,  we  need  to  determine  rn  such  that  |  Jny|  <  rB|y| 
for  vectors  y  in  S„,  and  tn  such  that  iJ^yl  ^  *»lvl  f°r  vectors  y  in  U„. 

Theorem.  (Sauer- Yorke,  1991)  Assume  6  <  ^-7,  where  m  is  the  dimension 
of  the  system,  and  let  B  be  a  bound  on  the  first  and  second  partial  derivatives 
of  the  map  f  and  its  inverse  /-1 .  Define  Co  =  0;  C„  =  esc  8  +  r„C„_i  for  n  >  0 
and  Dn  =  0;  Dn  =  csc#„  +  tnDn+i  for  n<  N.  If 

max{c"D")i^7s  <51) 

for  all  n  =  0, . . . ,  N,  then  there  exists  an  orbit  {y„}  of  f  such  that  |x„  —  y„|  < 
VS  forn  =  0,...,N. 

As  before,  8  is  the  noise  amplitude,  x„  is  the  noisy  orbit,  and  yn  is  the 
shadow  orbit.  The  time  step  N  is  considered  to  be  the  shadowing  time,  while 
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the  location  of  the  orbit  at  this  time  step  is  the  location  of  the  near-tangency. 
The  quantities  C„  and  Dn  are  obtained  through  recursion  relations  involving 
the  evolution  of  the  stable  and  unstable  manifolds  along  the  orbit  at  each  time 
step  n.  At  the  time-step  for  which  either  Cn  or  Dn  gets  larger  than  the  bound 
determined  by  the  theorem,  a  homoclinic  near-tangency  has  been  encountered 
and  the  system  is  no  longer  shadowable. 

The  nonhyperbolic  shadowing  theorem  contains  the  hyperbolic  shadowing 
theorem  of  Anosov  and  Bowen  as  a  specific  case;  for  hyperbolic  systems,  the 
quantities  Cn  and  Dn  remain  bounded  from  above,  and  never  violate  the  bound 
given  in  the  above  theorem.  Hence  for  hyperbolic  systems  the  shadowing  time 
N  is  infinite. 

The  Sauer  and  Yorke  theorem  also  quantifies  the  shadow  distance  e.  While 
in  the  Anosov- Bowen  theorem  the  shadow  distance  was  only  known  to  exist 
but  never  explicitly  specified,  in  the  new  theorem  it  is  seen  that  this  shadow 
distance  e  =  v6,  which  is  the  maximum  distance  between  the  noisy  and  deter- 
ministic orbits  as  long  as  the  orbit  is  not  integrated  beyond  n  =  N. 

5.4  The  Nonhyperbolic  Shadowing  Algorithm 
The  primary  method  for  finding  shadow  orbits  includes  the  refinement 
procedure  of  Hammel,  Yorke,  and  Grebogi  [38].  Given  a  noisy  and  deterministic 
trajectory,  the  idea  is  to  refine  the  noisy  orbit  by  applying  a  recursion  relation 
to  the  components  of  the  noisy  orbit  until  the  noise  can  no  longer  be  reduced. 
The  refinement  procedure  as  given  in  [38]  and  used  in  [39]  and  [40]  is  as 
follows. 
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Let  en+i  represent  the  one-step  error  for  the  noisy  orbit,  and  An  represent 
the  distance  between  the  noisy  orbit  and  a  deterministic  orbit  at  time  step  n: 

e„+i  =  y„+i  -  /(y»),  (5.2) 

A„  =  x„  -  y„.  (5.3) 

We  can  then  combine  these  equations  to  write 

An+i  =  /(x„)  -  /(y„)  -  en+i.  (5.4) 

We  now  assume  that  A„,  being  the  shadow  distance  for  a  given  time  step, 
is  small,  so  we  can  expand  /(xn )  in  a  Taylor  series  and  linearize  to  obtain 

/(x„)~/(y„)  +  J„An,  (5.5) 

where  Jn  is  the  Jacobian  of  the  appropriate  map;  if  the  dynamical  system  is  a 
system  of  differential  equations,  then  «/„  is  the  Jacobian  of  the  map  determined 
by  the  integration  routine.  So  equation  (5.4)  becomes 

An+i  =  J„A„  -  en+i.  (5.6) 

We  wish  to  find  {  A„}£_0  in  terms  of  the  coordinates  defining  the  stable  and 
unstable  manifold.  So  we  will  rewrite  An  and  e„  in  terms  of  the  variables  un 
and  sn.  These  are  the  unit  basis  vectors  which  lie  along  the  unstable  and  stable 
directions;  they  can  be  constructed  using  the  Gram-Schmidt  orthogonalization 
procedure.  So 

An  =  a„u„  +/?„s„,  (5.7) 

and 

en  =  *7i»u„  +  CnS„.  (5.8) 
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The  above  equations  can  then  be  substituted  into  equation  (5.6),  and  after 
some  algebraic  manipulation,  we  obtain  the  recursion  relations 

an  =  (an+i  +  »7n+l)/|  JnUn|,      "AT  =  0,  (5.9) 

and 

Pm+l=0u\J*»m\-G*+U     ^0  =  0.  (5.10) 

Computing  these  coefficients,  one  is  able  to  find  a  trajectory  that  is  less 
noisy  than  the  original  one.  Near  a  homoclinic  tangency,  however,  no  improve- 
ment over  the  noise  may  be  possible. 

Some  difficulties  are  inherent  in  this  procedure,  however.  First  of  all,  the 
refinement  computations  must  be  carried  out  at  a  higher  accuracy  than  the 
noise  level.  Thus,  if  it  is  desired  to  shadow  trajectories  with  double-precision 
level  computer  noise,  then  computations  using  96-bit  arithmetic  are  required. 
This  is  obviously  impractical.  An  added  complication  is  the  dependence  of  the 
shadow  distance  on  the  one-step  error.  It  thus  may  be  necessary  to  consider 
many  different  realizations  of  the  noise  before  a  reasonable  estimate  of  the 
shadow  distance  is  obtained. 

The  algorithm  obtained  from  the  new  shadowing  theorem  does  not  have 
these  difficulties.  First  of  all,  computations  need  only  be  carried  out  to  the 
same  precision  as  the  noise.  Also,  the  refinement  process  is  included  in  the 
theorem  itself,  where  the  bounds  on  the  noise  level  are  used  instead  of  the 
actual  one-step  error.  It  is  hence  not  necessary  to  know  the  exact  form  of  the 
noise  in  order  to  calculate  the  shadowing  time.  Of  course,  this  simplification  in 
the  numerical  process  of  calculating  shadowing  results  leads  to  the  limitation 
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of  only  obtaining  bounds  on  the  shadowing  time,  and  thus  the  results  obtained 
must  be  interpreted  somewhat  differently. 

Recall  from  chapter  two  that  invariant  manifolds  are  defined  anywhere 
within  a  hyperbolic  set.  Given  an  initial  condition  of  a  differential  equation,  the 
stable  and  unstable  manifolds  can  be  calculated  anywhere  along  the  trajectory, 
until  a  homoclinic  tangency  is  reached.  The  stable  and  unstable  directions, 
assuming  they  are  each  one- dimensional,  are  determined  by  the  equations 

U„+l  =  «7nU„/|J„Un|,  (5.11) 

Sn+1  =  /nS„/|J„Sn|.  (5-12) 

The  unstable  manifold  at  each  point  can  be  calculated  be  choosing  a  random 
vector  and  iterating  equation  (5.11)  forwards  in  time,  while  the  stable  manifold 
at  that  point  is  calculated  by  choosing  a  random  vector  and  iterating  equation 
(5.12)  backwards.  If  the  manifolds  are  more  than  one-dimensional,  than  Gram- 
Schmidt  orthogonalization  can  be  used. 


CHAPTER  6 

RESULTS  OF  SHADOWING 

NONHYPERBOLIC  SYSTEMS 

The  nonhyperbolic  shadowing  theorem  that  was  proven  by  Sauer  and  Yorke 
was  used  as  a  basis  for  numerical  experiments  involving  the  shadowability  of 
various  nonhyperbolic  systems.  These  numerical  experiments  involved  utilis- 
ing the  algorithm  outlined  in  the  previous  chapter,  primarily  to  calculate  the 
shadowing  time  and  the  location  in  the  system  where  shadowing  breaks  down. 
This  chapter  will  discuss  these  shadowing  results  for  a  number  of  nonhyperbolic 
systems. 

6.1  The  Model  Systems 
Shadowing  results  were  obtained  for  a  variety  of  systems  including  both 
maps  and  differential  equations.  The  maps  studied  are  the  Henon  map  and 
the  standard  map;  the  differential  equations  include  truncations  of  the  Toda 
lattice  potential  and  the  anisotropic  Kepler  potential.  While  the  potential 
systems  are  more  relevant  physically,  in  fact,  both  the  Toda  lattice  truncations 
and  the  anisotropic  Kepler  problem  have  applications  in  galactic  dynamics  and 
condensed  matter  physics,  the  maps  are  easier  to  analyze  and  computationally 
less  time-consuming.  Maps  are  also  relevant  for  study  because,  as  mentioned 
in  chapter  2,  they  can  be  considered  reasonable  approximations  to  Poincare 
maps  of  differential  equations. 
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Figure    6.1     The  Henon  map,  with  a=1.4  and  b=0.3. 


As  a  first  example,  consider  the  Henon  map    [41]: 


Zn+l  =hn+a-  xn 


J/n+l  =  Xn. 


(6.1) 


This  is  a  dissipative  map  as  long  as  b  is  chosen  to  be  less  than  one,  independent 
of  the  value  of  a.  In  this  study,  the  parameters  a  and  b  were  chosen  to  be  1.4 
and  0.3,  respectively,  so  that  iteration  of  the  map  yields  the  familiar  strange 
attractor  shown  in  figure  6.1. 

As  an  example  of  a  conservative  map,  the  shadowing  algorithm  was  applied 
to  the  standard  map,  which,  in  action-angle  variables,  has  the  form 


Jn+l  =  Jn  +  (K/2Tr)sin2ir0n  (mod  1), 


(6.2) 


#n+l  =  #n  +  Jn+1 


(mod  1). 


(6.3) 
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Figure    6.2     A  single  stochastic  trajectory  in  the  standard  map,  for  K  =  3. 


This  map  is  nonhyperbolic  for  all  positive  values  of  the  stochasticity  parameter 
K   [39]. 

The  standard  map  is  a  good  example  of  a  generic  nonhyperbolic  area- 
preserving  map,  conveniently  displaying  the  important  qualities  present  for 
such  systems.  For  a  typical  value  of  K,  the  phase  plane  is  divided  into  regions 
consisting  entirely  of  regular  orbits  and  regions  consisting  entirely  of  orbits 
that  wander  stochasticly  around  the  portion  of  the  phase  plane  accessible  to 
them.  These  regions  are  seperated  from  each  other  by  KAM  tori.  When  K  is 
small,  the  majority  of  the  phase  space  is  covered  by  regular  orbits,  with  isolated 
regions  of  stochasticity  associated  with  the  resonance  zones.  As  K  increases, 
the  measure  of  the  stochastic  region  of  phase  space  also  increases,  as  KAM 
tori  break  down  into  cantori,  and  the  cantori  themselves  begin  to  disappear. 
At  a  critical  value  of  K ,  appropriately  called  the  critical  value,  the  last  KAM 
torus  that  seperates  regions  of  stochasticity  disappears;  above  this  value  the 
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standard  map  exhibits  global  stochasticity.  A  typical  stochastic  orbit  is  shown 
in  figure  6.2  for  K  =  3,  which  is  above  the  critical  value.  The  regions  where 
the  orbit  does  not  visit  are  occupied  be  regular  orbits,  and  seperated  from  this 
stochastic  orbit  by  KAM  tori. 

The  first  differential  equation  studied  was  the  Hamiltonian  system  gener- 
ated by  the  Henon-Heiles  potential.  The  Henon-Heiles  potential  was  studied 
by  Henon  and  Heiles  [42]  in  the  1960's  as  a  simple  model  of  a  galaxy  which 
includes  stochastic  orbits.  This  is  a  two  degree-of-freedom  system  with  Hamil- 
tonian 

H  =  \(p2x+P2y)  +  V(x,y),  (6.4) 


with  the  potential  given  by 


1  ■  2  ,  „,2\  ,  Ji„      1    3 


V(»,  y)  =  ^(x2  +  y2)  +  x2y  -  -y\  (2.8) 

An  interesting  example  of  an  integrable  Hamiltonian  system  is  the  Toda 
lattice  [43].  The  periodic  three-particle  Toda  lattice  is  a  chain  of  three  particles 
connected  by  non- linear  springs,  and  is  given  by  the  Hamiltonian 

H  =  ~{qx2  +  <Z22  +  <Z32)  +  exp(9l  -  q2) 

&  (°-o) 

+  exp(g2  -  93)  +  exp(?3  -  Ql)  ~  3. 
Recall  that,  in  order  to  be  integrable,  this  system  needs  to  have  two  constants 
of  motion  besides  the  energy.  These  constants  have  been  proven  to  exist  and 
have  been  found  by  Flaschka  [44]  and  Henon  [45].  Flaschka's  method  involved 
the  use  of  Lax  pairs  [46].  A  Lax  pair  is  a  pair  of  matrices,  say  L  and  M,  whose 
elements  depend  on  a  parameter  t  such  that  ^  =  LM  -  ML.  The  constants  of 
the  motion  are  then  the  coefficients  in  the  charactristic  polynomial  det(L  —  XL), 
which  do  not  vary  with  t. 
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The  three-particle  Toda  lattice  can  be  reduced  to  a  two-degree  of  freedom 
system  via  a  transformation 


(6.6) 


H  =  l-{i2  +  y2)  +  ^[exp(2y/3x  +  2y) 

+  exp(-2VZx  +  2y)  +  exp(-4y)]  -  -. 

Expanding  the  exponential  terms  up  to  sixth  order,  we  obtain 
tt        1 <  -2    ,     -2    ,      2    ,      2\    ,      2  *    3 

=  2~(x  +  v  + x  +  y)  +  *»  "  ^ 

+  ^4  +  x2y2  +  \yA  +  x4y  +  ^x2y3  -  \y5  (6.7) 

The  second-order  truncation  of  this  potential  in  integrable,  consisting  sim- 
ply of  two  linear  harmonic  oscillators.  The  third-order  truncation  is  seen  to 
be  the  Henon-Heiles  potential,  which  is  nonintegrable.  In  fact,  it  has  been 
proven  [47]  that  all  truncations  of  the  periodic  three-particle  Toda  lattice 
above  second-order  are  nonintegrable  and  include  both  regular  and  stochastic 
orbits.  The  specific  case  of  the  sixth-order  truncation  of  this  Toda  lattice  will 
also  be  studied. 

The  truncations  of  the  Toda  lattice  present  relatively  similar  shadowing 
problems  in  that  the  potentials  only  involve  polynomials.  In  contrast,  shadow- 
ing of  the  anisotropic  Kepler  potential  will  also  be  attempted.  In  normalized 
cartesian  coordinates,  the  Hamiltonian  for  the  anisotropic  Kepler  potential  in 
two  dimensions  is  given  by 

22  -i 

^  =  ^  +  J___J=,  (6.8) 

which  corresponds  to  a  system  with  no  angular  momentum,  and  where  u  and 
v  are  the  momenta.  Besides  the  energy,  this  system  is  also  characterized  by 
another  free  parameter,  the  mass  ratio  fi/u. 
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Shadowing  of  the  anisotropic  Kepler  problem  is  more  difficult  than  the  pre- 
vious problems  studied  primarily  due  to  the  singularity  at  the  origin.  Because 
of  this,  the  bounds  required  to  calculate  the  quantities  in  the  nonhyperbolic 
shadowing  theorem  are  no  longer  defined;  hence,  a  new  modification  of  the 
theorem  is  required. 

6.2  Shadowing  Time  Results  for  Maps 

For  all  of  the  systems  discussed  in  section  6.1,  the  shadowing  time  was 
calculated  given  the  noise  amplitude  for  a  large  number  of  initial  conditions. 
These  numerical  calculations  involved  computing  the  quantities  associated  with 
the  nonhyperbolic  shadowing  theorem  discussed  in  chapter  5  until  the  values 
for  C„  and  Dn  exceeded  their  bounds,  at  which  point  the  shadowing  theorem 
breaks  down  and  the  trajectory  is  no  longer  shadowable.  The  locations  in  phase 
space  where  this  happens  are  the  homoclinic  near-tangencies  of  the  system, 
hereafter  called  'glitches.'  Recall  from  chapter  5  that  the  shadowing  time  is  the 
length  of  time  that  the  trajectory  is  computed  before  a  glitch  is  encountered. 
For  each  noise  amplitude  8,  the  location  of  the  glitch  and  the  shadowing  time 
were  determined  for  an  ensemble  of  1000  randomly  chosen  initial  conditions. 
For  all  systems,  noise  amplitudes  were  chosen  to  be  from  8  =  10  to  8  =  10  . 
Since  all  important  quantities  were  calculated  in  double  precision,  8  =  10-14, 
which  roughly  corresponds  to  double  precision  error,  is  the  lower  limit  on  the 
noise  amplitude  computationally. 

The  first  system  studied  was  the  Henon  map.  The  shadowing  time  was 
computed  for  each  of  the  1000  randomly  chosen  initial  conditions  for  each  noise 
level.  Figure  6.3  shows  the  distribution  of  shadowing  times  for  representative 
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Figure  6.3  Histograms  showing  the  distribution  of  shadowing  times  for 
randomly  chosen  ensembles  of  1000  initial  conditions,  chosen  to  lie  within  10-5 
of  the  origin,  for  the  Henon  map.  (a)  8  =  10-14;  (b)  6  =  10-13; 
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Figure  6.3.        (continued) 

(c)£=10-12;(d)<S  =  10-n; 
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Figure  6.3.        (continued) 
(e)  6  =  10-10;  (f)  6  =  10"8. 
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ensembles  of  1000  initial  conditions,  chosen  to  lie  within  a  distance  of  10-5  to 
the  origin. 

The  first  conclusion  that  can  be  inferred  from  this  data  is  that  the  shadow- 
ing time  varies  significantly  for  this  system,  spanning  several  orders  of  magni- 
tude. The  conventional  wisdom  is  that  the  shadowing  time  for  typical  nonhy- 
perbolic  systems  is  long  enough  to  consider  any  deviation  from  hyperbolicity 
to  be  insignificant.  However,  we  can  see  that  even  though  the  mean  shadowing 
time  can  be  long,  there  is  an  appreciable  number  of  orbits  with  relatively  short 
shadowing  times,  extending  down  to  about  two  orders  of  magnitude  below  the 
mean. 

Second,  these  distributions  remain  essentially  identical  if  one  chooses  a 
different  ensemble  of  initial  conditions,  for  example,  by  grouping  the  initial 
conditions  around  a  randomly  chosen  location  instead  of  the  origin,  or  changing 
the  range  around  the  location  about  which  the  initial  conditions  are  chosen  to 
he. 

This  last  observation  is  of  extreme  importance;  it  seems  to  imply  that 
the  dependence  of  shadowing  time  on  initial  conditions  is  fractal.  In  other 
words,  the  distributions  have  similar  structure  on  all  scales.  Choosing  the 
initial  conditions  randomly  within  a  region  10~10  from  the  origin  yields  sta- 
tistically near-identical  distributions  as  choosing  the  initial  conditions  within 
a  region  10-5  from  the  origin.  This  conclusion  is  supported  by  results  of  the 
Kolmogorov-Smirnoff  test,  which  shows  that  these  sample  distributions  have 
about  a  95  percent  chance  of  being  drawn  from  the  same  distribution. 

By  comparing  the  sample  distributions  in  figure  6.3,  it  is  apparent  that  the 
mean  shadowing  time  decreases  as  the  noise  level  is  increased.  This  can  easily 
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Figure  6.4  log  ns  versus  log  6  for  the  Henon  map,  where  ns  is  the  mean 
shadowing  time  and  8  is  the  noise  level.  The  error  bars,  obtained  by  averaging 
over  five  distributions,  are  smaller  than  the  squares. 
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be  understood  -  increasing  the  noise  level  increases  the  range  of  angles  which 
are  close  enough  to  zero  to  create  a  homoclinic  near-tangency,  making  the  orbit 
experience  these  glitches  more  frequently.  The  relationship  between  the  mean 
shadowing  time,  ns,  and  the  noise  level  8  is  shown  in  figure  6.4.  Each  point 
was  determined  by  averaging  over  five  distributions.  Since  the  distributions  are 
so  close  to  being  identical  statistically,  the  standard  deviation  of  each  point  is 
smaller  than  the  squares  used  for  plotting. 

Hammel,  Yorke  and  Grebogi  predict  a  power  law  relationship  between  the 
shadowing  time  and  the  noise  amplitude,  as  provided  in  the  following  conjec- 
ture   [38]: 

Conjecture.  For  a  typical  dissipative  map  f  :  R2  — ►  R2  with  a  positive  Lya- 
punov  exponent  and  a  small  noise  amplitude  8  >  0,  we  expect  to  End  e  <  v8 
for  an  orbit  length  N  ~  \/\/8. 

This  conjecture  is  not  very  meaningful  when  applied  to  individual  orbits, 
due  to  the  fractal  distributions  of  shadowing  times.  However,  if  one  interprets 
the  conjecture  instead  as  applying  to  the  mean  shadowing  time,  then  it  implies 
effectively  that  ns  ~  1/V8  for  dissipative  maps. 

As  seen  in  figure  6.4,  The  six  data  points  for  8  =  10-14  through  8  =  10-9 
are  well  fit  by  a  power  law  ns  ~  6P,  where  p  =  —0.556  ±  0.006.  Since  the 
standard  deviation  for  each  point  is  insignificant,  the  error  estimate  for  the 
exponent  p  is  simply  the  formal  error  associated  with  the  least  squares  fit  of 
the  average  values  plotted  in  figure  6.4. 

These  results  roughly  support  the  scaling  conjectured  by  Hammel,  Yorke 
and  Grebogi,  and  a  very  precise  scaling  law  is  evident,  though  it  is  not  exactly 
equal  to  the  law  n3  ~   l/v8  predicted,  at  least  within  the  range  of  noise 
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amplitudes  8  =  10-14  through  8  =  10-9.  If  one  were  to  include  the  data  points 
for  8  =  10-8  and  8  =  10-7,  the  exponent  would  be  p  =  -0.624  ±  0.030.  This 
trend  towards  slightly  shorter  shadowing  times  than  expected  for  high  noise 
amplitude  appears  to  be  due  to  the  'spikes'  in  the  histograms  for  these  noise 
amplitudes. 

The  fractal  nature  of  the  shadowing  time  distributions  appear  to  be  a 
universal  feature,  which  is  not  surprising.  Similar  distributions  are  observed 
in  the  standard  map,  and  are  shown  in  figure  6.5  for  the  case  that  K  =  7.  At 
large  values  of  K,  the  measure  of  the  regular  region  goes  to  zero.  Thus,  this 
value  of  K  was  chosen  to  ensure  that  all  orbits  would  be  within  the  stochastic 
sea,  and  that  there  would  be  no  effect  due  to  choosing  initial  conditions  close 
to  regular  regions  of  phase  space.  Again,  these  distributions  are  independent  of 
the  range  and  location  of  the  chosen  initial  conditions.  These  results  are  gain 
supported  by  the  Kolmogorov- Smirnoff  test,  which  verifies  that  these  sample 
distributions  have  a  probability  of  within  90%  of  being  drawn  from  the  same 
distribution.  Similar  results  are  obtained  for  the  value  K  =  3. 

A  conjecture  similar  to  that  for  dissipative  maps  has  been  reported  by 
Grebogi  et  al    [39]  for  Hamiltonian  maps: 

Conjecture.  For  a  typical  two-dimensional  Hamiltonian  map  yielding  chaotic 
trajectories  with  a  small  noise  amplitude  8  >  0,  we  expect  to  find  e  <  V8  for 
a  trajectory  of  length  N  ~  1/Vti. 

The  relationship  between  the  mean  shadowing  time,  ns,  and  the  noise  level 
8  is  shown  in  figure  6.6.  Again,  the  data  points  were  obtained  by  averaging  ns 
over  five  different  ensembles  of  1000  initial  conditions  each  for  each  noise  level. 
Fitting  the  data  points  to  a  straight  line  yields  ns  ~  8P,  where  p  =  -0.498  ± 
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Figure  6.5  Histograms  showing  the  distribution  of  shadowing  times  for 
randomly  chosen  ensembles  of  1000  initial  conditions,  chosen  to  lie  within  10-5 
of  the  origin,  for  the  standard  map  with  K  =  7.  (a)  6  =  10~14;  (b)  6  =  10~12; 
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(c)  6  =  10"10;  (d)  6  =  10"8. 
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0.002,  in  complete  agreement  with  the  conjecture  given  the  interpretation  that 
the  trajectory  length  should  be  the  mean  shadowing  time. 

6.3  Shadowing  Time  Results  for  Hamiltonian  Differential  Equations 
Similar  shadowing  time  distributions  are  observed  for  the  case  of  chaotic 
Hamiltonian  differential  equations,  as  observed  for  the  Henon-Heiles  potential 
(figure  6.7)  given  an  energy  E  =  0.16  and  for  the  sixth-order  truncation  of  the 
Toda  Lattice  (figure  6.8)  for  an  energy  of  E  =  30.  In  each  case,  the  energy 
was  chosen  to  yield  a  significant  measure  of  stochastic  orbits.  For  the  Henon- 
Heiles  potential,  results  were  obtained  for  an  energy  of  E  =  0.14  as  well,  with 
statistically  similar  results. 

The  shadowing  method  for  two  degree-of-freedom  Hamiltonian  systems  is 
similar  to  the  method  used  for  two-dimensional  maps.  The  primary  difficulty 
is  due  to  the  calculation  of  the  stable  and  unstable  manifolds.  A  significantly 
larger  number  of  time  steps  is  required  for  the  manifolds  to  converge  for  the  case 
of  differential  equations  and,  in  fact,  they  are  often  observed  not  to  converge  at 
all,  in  a  reasonable  amount  of  time.  This  is  indeed  the  case  for  the  sixth-order 
truncation  of  the  Toda  lattice;  if  the  energy  is  increased  much  beyond  E  =  30, 
the  phase  space  becomes  more  and  more  stochastic,  and  the  stable  manifold  is 
not  easily  computed.  Even  for  E  =  30,  many  choices  of  initial  conditions  did 
not  yield  reliable  results  due  to  the  inability  to  calculate  the  stable  manifold, 
and  so  the  distributions  shown  in  figure  6.8  are  for  a  smaller  number  of  initial 
conditions. 

While  generic  two  degree-of-freedom  systems  have  a  four- dimensional  phase 
space,  two  degree-of-freedom  Hamiltonian  systems  have  only  one- dimensional 
stable  and  unstable  manifolds.  Due  to  Liouville's  theorem,  the  number  of  stable 
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Figure  6.6  log  n3  versus  log  6  for  the  standard  map  with  K  =  7.  The 
error  bars,  obtained  by  averaging  over  five  distributions,  are  smaller  than  the 
squares. 
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(contracting)  directions  must  equal  the  number  of  unstable  (expanding)  direc- 
tions, to  preserve  phase  space  volume.  Since  E  =  const,  can  be  used  to  reduce 
the  phase  space  dimension  by  one,  and  the  direction  along  the  trajectory  can 
be  neither  expanding  nor  contracting,  the  stable  and  unstable  manifolds  can 
each  have  only  one  dimension.  The  numerical  method  of  calculating  shadowing 
times  in  this  case  is  shown  in  appendix  A.  If  the  system  is  not  Hamiltonian, 
these  manifolds  are  not  necessarily  one-dimensional  and,  in  fact,  the  dimension 
of  the  manifolds  may  actually  change  with  time.  This  case  is  considered  by 
Dawson  et  al  [48],  who  determine  that  the  shadowing  time  can  actually  be 
much  lower  in  this  case,  due  to  the  fluctuating  number  of  dimensions  of  the 
invariant  manifolds.  The  calculation  of  shadowing  times  for  two-dimensional 
stable  and  unstable  manifolds  can  be  done  using  the  Fortran  code  shown  in 
appendix  B. 

Unlike  the  case  for  maps,  the  relation  between  the  mean  shadowing  time 
and  the  noise  amplitude  is  not  a  simple  power  law  for  these  potentials.  This 
is  shown  in  figure  6.9  and  figure  6.10.  This  should  not  be  too  surprising,  since 
there  is  no  a  priori  reason  why  this  relation  should  be  a  power  law.  Recall  that 
a  'glitch'  occurs  when  the  angle  between  the  stable  and  unstable  manifolds  is 
small  enough  that  the  noise  amplitude  is  able  to  cause  a  trajectory  to  'jump 
over'  a  stable  manifold.  It  is  therefore  apparent  that  the  mean  shadowing  time 
should  be  a  monotonically  decreasing  function  of  the  noise  amplitude.  There 
is  no  guarantee,  however,  that  this  relationship  be  a  simple  function,  e.g.  a 
power  law.  Indeed,  as  a  particle  moves  along  a  trajectory,  the  angle  between 
the  stable  and  unstable  manifolds  fluctuates  smoothly  in  the  case  of  continuous 
systems,  but  the  angle  is  not  necessarily  a  monotonic  function  of  time.  As  this 


90 


100 


T 1     I    I   I  I  III 1 1 — I    I   I  I  II 


I     I    I   I  I  III 1 1     I    I  I  I  II 


90  [  (a)  6  =  10 

80 
70 
60 
50 
40 
30 
20 
10 
0 


-14 


Si 


1 


J I I  I  I  I  III I I I  I  I  I II 


I Il I I I  I  I  I II 


.01 


.1  1  10 

Shadowing  Time,  in  crossing  times 


.01 


100 


uu 

1      1     1    1  1  1  ll|           1       1     1    I  I  I  1 1 1           1       1     1    1  1  1  llj           1 

1    1  1  II II 

90 

— 

(b)  6  =   10"13 

— 

80 

— 

70 

— 

— 

60 

— 

50 

n 

— 

40 
30 

tJ 

w 

\ 
ftf 

J 

— 

20 

— 

— ' 

— 

10 

{ 

0 

I      i    I 

i  i  i  ill          i 

I    I 

i 

ill 

1 

1    1 

1  1  1 1 1           1 

1    1  1  1  1 1 1 

.1  1  10 

Shadowing  Time,  in  crossing  times 

Figure  6.7  Histograms  showing  the  distribution  of  shadowing  times  for 
randomly  chosen  ensembles  of  1000  initial  conditions,  chosen  to  lie  within  10-5 
of  the  origin,  for  the  Henon-Heiles  potential  with  E  =  0.16.  (a)  6  =  10~14;  (b) 
6  =  10"13; 
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Figure  6.7.        (continued) 

(c)^  =  10-n;(d)^=10-10. 
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angle  fluctuates,  it  then  has  several  maxima  and  minima;  this  causes  certain 
noise  amplitudes  to  intersect  the  angle  versus  time  curve  more  frequently  that 
others. 

The  power  law  conjectures  of  Grebogi,  Hammel,  Sauer  and  Yorke  actually 
hold  up  quite  well  in  the  case  of  these  continuous  systems  if  only  a  smaller  range 
of  noise  amplitudes  is  considered,  for  example,  from  8  =  10-12  to  8  =  10-10  in 
these  cases.  For  very  low  noise  amplitudes,  however,  the  shadowing  times  are 
considerably  shorter  than  one  would  predict  from  the  scaling  law  ns  ~  8P. 

6.4  Results  Concerning  Glitch  Locations 
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Figure  6.11     The  glitch  locations  for  the  Henon  map  for  an  ensemble  of  1000 
initial  conditions  and  8  =  10-14. 
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Figure  6.8  Histograms  showing  the  distribution  of  shadowing  times  for 
randomly  chosen  ensembles  of  1000  initial  conditions,  chosen  to  lie  within  10-5 
of  the  origin,  for  the  sixth-order  Toda  potential  with  E  =  30.  (a)  6  =  10~14; 
(b)  6  =  lO"13; 
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Figure  6.8.        (continued) 

(c)£  =  10-12;(d)£=10-n. 
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Figure   6.9     log  ns  versus  log  6  for  the  Henon-Heiles  potential  with  E  =  0.16. 
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Figure    6.10     log  n3  versus  log  6  for  the  sixth-order  truncation  of  the  Toda 
lattice  with  E  =  30. 
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Figure    6.12     The  glitch  locations  for  the  standard  map  with  K  =  3  for  an 
ensemble  of  1000  initial  conditions  and  6  =  10-14. 


Recall  that  a  glitch  is  defined  to  be  the  location  at  which  the  stable  and 
unstable  manifolds  are  close  enough  to  being  tangent  that  the  trajectory  is 
not  shadowable  at  that  point  for  the  given  noise  amplitude.  For  each  initial 
condition,  the  glitch  locations  were  determined  as  well  as  the  shadowing  time. 
For  the  Henon  map,  the  locations  of  the  glitches  for  one  representative  ensemble 
of  initial  conditions  with  6  =  10-14  are  shown  in  figure  6.11.  It  should  be 
apparent  that,  if  one  were  to  continue  computing  the  locations  of  these  near- 
tangencies,  the  resulting  figure  would  not  be  distinguishable  from  the  Henon 
attractor  itself,  figure  6.1.  Thus,  it  appears  that  the  homoclinic  near-tangencies 
are  dense  within  the  attractor.  Similar  results  are  found  for  other  choices  of 
the  noise  level. 
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Figure    6.13     The  glitch  locations  for  the  standard  map  with  K  =  7  for  an 
ensemble  of  1000  initial  conditions  and  8  =  10-14. 


The  situation  is  somewhat  different  for  area-preserving  maps,  which  do  not 
have  attractors.  In  figure  6.12,  the  glitch  locations  are  shown  for  the  standard 
map  with  K  =  3.  As  can  be  seen,  these  glitches  are  not  randomly  distributed 
throughout  the  stochastic  region  of  phase  space,  but  appear  to  lie  along  distinct 
curves  within  the  stochastic  region.  Comparing  the  glitch  locations  to  the 
stochastic  region  of  phase  space  observed  in  figure  6.2,  it  is  observed  that  a 
number  of  these  glitches  occur  at  the  boundary  of  the  stochastic  and  regular 
regions. 

A  similar  distribution  of  glitches  occurs  for  K  =  7  as  well,  as  shown  in 
figure  6.13,  even  though,  as  shown  in  figure  6.14,  practically  the  entire  phase 
space  appears  stochastic. 
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Figure    6.14     A  single  stochastic  orbit  for  the  standard  map  with  K  =  7. 


CHAPTER  7 
DISCUSSION  AND  CONCLUSIONS 

This  chapter  is  devoted  to  a  recapitulation  of  the  major  results  of  this 
study,  as  well  as  some  general  remarks  concerning  the  interpretation  of  these 
results. 

The  results  of  this  study  seem  to  indicate  that,  even  for  gravitational  sys- 
tems with  a  large  number  of  particles,  discrete  effects  due  to  close  encounters 
between  particles  have  an  important  effect  on  the  instability  of  trajectories  in 
the  system.  These  results  hold  provided  that  the  system  of  interest  is  nonin- 
tegrable  and  allows  for  a  finite  fraction  of  stochastic  orbits.  Numerical  studies 
seem  to  indicate  that  this  is  in  fact  the  case  for  elliptical  galaxies. 

These  discreteness  effects  are  seen  to  increase  the  rate  of  instability  in  grav- 
itational TV-body  systems,  as  discussed  in  chapter  3.  Results  of  these  numerical 
studies  indicate  that  changes  in  the  softening  parameter,  which  was  used  to 
alter  the  effect  of  close  encounters,  certainly  changed  the  rate  of  instability. 
There  is  no  sense  in  which  the  instability  "turns  off'  when  N  is  large,  and  close 
encounters  remain  important.  Our  numerical  studies  hence  indicate  that,  even 
for  a  large  number  of  stars,  both  the  effects  of  the  smooth  mean  field  and  the 
close  encounter  noise  are  important  when  determining  the  rate  of  instability 
in  gravitational  N-body  systems. 

Analytical  studies  also  indicate  how  a  nonintegrable  mean  field  and  noise 
due  to  discreteness  effects  are  both  important  when  determining  the  rate  of 


100 


101 
instability.  The  primary  conclusion  of  chapter  4  is  that  the  presence  of  a  nonin- 
tegrable  mean  field  drastically  alters  the  relaxation  time,  which  was  calculated 
by  Chandrasekhar  in  the  idealised  case  that,  without  collisions,  the  stars  fol- 
low a  smooth,  rectilinear  trajectory.  These  new  calculations  indicate  that  the 
exponential  instability  of  stellar  trajectories  in  the  nonintegrable  mean  field 
reduces  the  relaxation  time  significantly. 

The  nonhyperbolic  shadowing  theorem  was  used  to  determine  how  large 
an  effect  noise  of  a  given  amplitude  may  have,  given  that,  in  most  physically 
relevant  situations,  the  nonintegrable  mean  field  is  nonhyperbolic.  Results 
from  both  maps  and  differential  equations  support  the  claim  that  systems 
which  include  dynamical  noise  can  be  shadowed  for  long  times  for  most  orbits, 
provided  that  the  noise  amplitude  8  is  small  enough.  However,  the  existence 
of  orbits  with  appreciably  shorter  shadowing  times  should  cause  concern  when 
integrating  chaotic  dynamical  systems  for  long  times.  The  existence  of  these 
'unshadowable'  orbits  makes  conclusions  about  the  'shadowability'  of  orbits  in 
a  nonhyperbolic  system  unclear,  even  though  the  mean  shadowing  time  may  be 
quite  large.  This  is  especially  true  considering  the  fractal  nature  of  the  sample 
distributions  of  shadowing  times.  If  one  is  interested  in  the  shadowing  time 
for  an  orbit  with  a  certain  initial  condition,  one  has  to  realize  that  there  are 
orbits  with  initial  conditions  arbitrarily  close  to  the  desired  initial  condition 
with  drastically  different  shadowing  times. 

However,  nonhyperbolic  shadowing  may  be  used  to  determine  the  fraction 
of  orbits  that  may  be  qualitatively  changed  by  the  presence  of  noise.  Even 
given  the  fractal  nature  of  the  shadowing  times,  the  sample  distributions  yield 
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information  as  to  how  many  orbits  are  'unshadowable'  given  the  number  of 
iterations  of  the  map. 

As  an  example,  consider  the  situation  in  which  one  is  concerned  with  how 
computer  errors  are  going  to  affect  orbits  in  a  deterministic  system.  We  can 
see  from  the  histograms  in  chapter  6  that,  even  for  double  precision  calcu- 
lations, the  majority  of  orbits  in  the  Henon-Heiles  potential  and  sixth-order 
Toda  lattice  truncation  are  unshadowable  after  only  ten  crossing  times,  and 
several  are  unshadowable  after  less  than  one  crossing  time.  In  the  case  of  the 
Henon-Heiles  potential,  the  fraction  of  orbits  that  become  unshadowable  in 
less  than  a  crossing  time  is  about  ten  percent.  Thus,  any  orbit  that  behaves 
peculiarly,  even  within  the  first  crossing  time,  should  be  considered  suspect. 
Since  these  distributions  appear  fractal,  it  is  impossible  to  determine  if  such 
unusual  orbits  are  actual  orbits  of  the  deterministic  system,  or  are  associated 
with  a  breakdown  in  shadowing,  even  in  principle. 

It  remains  to  be  seen  whether  'unshadowability'  of  orbits  actually  affects 
the  statistical  properties  of  the  orbits  involved,  however.  It  is  certainly  true 
that  if  a  system  is  shadowable,  then  it  is  structurally  stable  with  respect  to 
noise.  In  fact,  hyperbolic  systems  have  been  proven  to  be  structurally  stable 
by  Anosov  [49].  The  converse  of  this  statement,  however,  is  not  necessarily 
true.  In  fact,  for  the  standard  map  with  large  stochasticity  (K  =  7),  the 
existence  of  homoclinic  near-tangencies  does  not  appear  to  alter  the  statistical 
properties  of  the  orbits.  This  would  seem  to  imply  that  the  standard  map 
for  K  =  7  is  structurally  stable  in  the  sense  that  the  statistical  properties  of 
the  map  is  unaffected  by  the  presence  of  noise.  It  remains  to  be  seen  whether 
this  behavior  is  universal.  For  example,  it  has  been  observed  that  even  weak 
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noise  can  affect  the  near-invariant  measure  of  Hamiltonian  potentials     [50]. 
One  would  expect  that  this  behavior  is  due  to  the  existence  of  homoclinic 
near-tangencies  within  the  system. 

The  fact  that  the  glitches  in  the  standard  map  are  not  randomly  distributed 
would  seem  to  imply  that  the  location  of  the  homoclinic  near-tangencies  may 
be  related  to  some  other  property  of  the  dynamical  system,  for  example,  the 
Lyapunov  exponents.  The  existence  of  glitches  adjacent  to  KAM  tori  may 
be  associated  with  'sticky'  orbits  with  smaller  local  Lyapunov  exponents  than 
are  found  in  the  stochastic  sea.  These  sticky  orbits  have  been  observed  for  the 
standard  map  [51]  and  for  several  two-dimensional  potentials  as  well  [52].  This 
may  have  important  consequences;  orbits  which  experience  a  glitch  adjacent 
to  invariant  tori  or  cantori  may  perhaps  'jump'  across  the  barrier,  even  if  a 
purely  deterministic  orbit  would  not.  These  glitches  would  then  serve  as  the 
mechanism  for  diffusion  across  these  barriers  in  noisy  dynamical  systems. 

It  is  important  to  note  that  the  results  concerning  the  shadowing  times  do 
not  seem  to  depend  on  the  form  of  noise  present.  It  is  certainly  true  that  for 
hyperbolic  systems  the  form  of  noise  does  not  really  matter  when  determining 
the  shadowability  of  the  system;  the  form  of  noise  will  alter  the  details  of 
individual  orbits,  but  the  statistical  properties  of  a  hyperbolic  system  will 
remain  unaltered,  since  they  are  structurally  stable  in  the  presence  of  noise. 
For  nonhyperbolic  systems,  however,  the  situation  is  more  complicated. 

The  shadowing  results  mentioned  in  the  previous  chapter  appear  to  be 
independent  of  the  form  of  noise;  only  the  amplitude  of  the  noise  was  considered 
to  be  a  factor.  The  important  thing  to  realize  is  that  the  shadowing  time  must 
be  interpreted  as  an  upper  bound  on  the  effect  the  noise  will  have;  A  noisy  orbit 
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is  not  guaranteed  to  be  affected  at  the  homoclinic  near-tangency  it  encounters; 
it  is  only  guaranteed  not  to  be  affected  if  it  does  not  approach  a  near-tangency. 
A  certain  form  of  noise  may  have  less  of  an  effect  than  this  bound  indicates; 
however,  as  long  as  the  amplitude  of  the  noise  is  given,  the  shadowing  theorem 
guarantees  that  the  choice  of  a  certain  form  of  noise  will  not  have  more  of  an 
effect  that  indicated. 

It  could  also  be  argued  that,  even  if  an  orbit  is  unshadowable,  the  statistical 
properties  of  the  system  may  remain  unaffected,  as  it  appears  for  the  standard 
map.  Assume  a  noisy  orbit  encounters  a  homoclinic  near-tangency  at  time  step 
JV,  after  which  the  orbit  no  longer  is  shadowed  by  the  original  deterministic 
orbit.  This  orbit  will  then,  in  all  probability,  begin  to  be  shadowed  by  a  new 
deterministic  orbit  beginning  at  time  step  TV  +  1.  Despite  the  fact  that  this 
noisy  orbit  will  not  be  similar  to  any  purely  deterministic  orbit  of  the  system, 
statistical  tests  will  most  likely  indicate  no  difference  between  the  noisy  and 
deterministic  orbits. 

The  only  situation  for  which  a  noisy  orbit  will  not  yield  qualitatively  similar 
results  to  what  is  expected  is  when  the  orbit  passes  between  a  stochastic  region 
of  phase  space  and  a  regular  region,  i.e.  that  it  crosses  KAM  tori  in  the  process. 
In  this  case,  the  qualitative  behavior  of  the  orbit  is  affected.  It  may  be  expected 
that  this  should  not  happen  very  often  for  systems  with  global  stochasticity,  in 
which  most  orbits  are  well  within  the  stochastic  sea,  far  away  from  any  regular 
regions  remaining.  However,  as  observed  in  the  case  of  the  standard  map,  the 
location  of  homoclinic  near-tangencies  seem  to  be  correlated  with  the  presence 
of  KAM  tori,  or,  at  least,  cantori.  Thus,  the  effect  of  these  near-tangencies 
may  be  more  pronounced  that  how  it  originally  appears. 


APPENDIX  A 
Fortran  Codes  -  two-dimensional  systems 

Program  Shadow 

This  program  calculates  the  shadowing  tine  and  location  of  glitch 
for  a  two  degree-of -freedom  Hamiltonian  differential  equation, 
specifically,  the  sixth-order  truncation  of  the  Toda  lattice. 
1  list  of  1000  initial  conditions  is  randomly  determined  and 
stored  in  a  file.    The  shadowing  times  for  each  initial  condition 
is  then  stored  in  a  file,  and  the  location  of  the  glitch  for  each 
initial  condition  is  stored  in  a  third  file.    The  map,  its 
Jacobian,  and  the  inverse  of  the  Jacobian  must  be  input  as 
seperate  subroutines. 

Variables  — 

xl,yl,xdl,ydl  —  the  initial  conditions 

x2,y2,xd2,yd2  —  the  glitch  location 

ran2  —  random  number  function,  from  Numerical  Recipes 

idum  —  random  number  seed 

i  —  index  to  calculate  ith  initial  condition 

ltime  —  the  shadowing  time 

lreg  —  Warning  flag  for  regular  or  escaping  point 

Files  ~ 

filel  —  the  initial  conditions 

file2  —  the  initial  conditions  and  the  shadowing  time,  or 

a  warning  statement 
file3  —  the  glitch  location 

todadat  —  the  three  new  files,  followed  by  the  model 
parameters 

Double  precision  xl,yl,xdl,ydl,x2,y2,xd2,yd2 

double  precision  sigma 

double  precision  delta, E.diaen 

Real  ran2 

Integer  i,j ran, lreg, ltime 

data  idum/-3303039/ 

character*12  filel, file2,file3 

common  //  delta, E.dimen 

open (unit=25 , f  ile= ' todadat ' , status= ' old ' ) 

read(25,*)filel,file2,file3 

open  (unit2 18, f ile=filel,status='new') 

open  (unit=19,f ile=file2,status='new') 

open  (unit =20, f ile=file3,status='new') 

open  (unit=21 ,f ile='smf ld.dat ' ,status='old') 

open  (unit=22,f ile= 'umfld.dat' ,status='old') 

open  (unit=23,f ile=' coords.dat' , status =' old') 
* 
*     Initialize  and  write  initial  conditions. 

read(25,*)delta,E,dimen 
do  i=0,999 
xl=0.d0 
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yl=ran2(idum) 
ydl=ran2  (  iduin) 

xdl=dsqrt (2 . dO*E-ydl*ydl-yl*yl+2 . dO*y l*yl*yl/3 . dO 
!  -yl**4+2.d0*yl**5/3.d0-22.d0*yl**6/45.d0) 

write(18,*)xl,yl,xdl,ydl 
end  do 

rewind  (unit=18) 
* 

*  For  each  initial  condition,  call  the  main  subroutine,  and  store 

*  results. 
* 

do  ii=l,1000 

ltime=0 

lreg=0 

read(18,*)xl,yl,xdl,ydl 

call  mfld(xl,yl,xdl,ydl,x2,y2,xd2,yd2,ltime,lreg,sigma) 

if  (lreg  . eq.   l)then 
ltime=0 

endif 

write(19,*)ii,xl,ydl,ltime,lreg 

write(20,*)x2,y2,xd2,yd2 
end  do 

close (unit = 18 , status= 'delete ' ) 
end 

Subroutine  mfld(xl,yl,xdl,ydl,x2,y2,xd2,yd2,ltime,lreg,sigma) 

Given  the  initial  conditions  xl,yl,xdl,ydl,  compute  the  shadowing 
time  and  the  location  of  the  glitch  for  the  sixth-order  trun- 
cation of  the  Toda  lattice.    Return  these  values  or  an  error  flag 
to  the  driver  program.    The  noise  level  delta  and  the  energy  must 
be  specified. 

Variables  — 

x(n),y(n),xdl(n),ydl(n)  —  The  phase  space  variables  at  time 

step  n 
z  —  swap  storage  used  when  computing  x(n+l) ,y(n+l) 
ajac(i.j)  —  the  jacobian  elements  of  the  standard  map 
ajacinv(i.j)  —  the  inverse  jacobian  elements 
anorm(i)  —  the  vector  used  to  normalize 
u(i,n) ,s(i,n)  —  the  stable  and  unstable  manifolds,  ith 

component,  at  time  step  n 
csc(n)  —  the  cosecant  of  the  angle  at  time  step  n 
csec  —  the  cosecant  for  the  current  time  step 
ul,u2,sl,s2  —  the  two  elements  of  the  stable  and  unstable 

manifolds  at  the  current  time  step 
m(n),tn(n)  —  the  bounds  on  the  expansion  of  the  map  at 

time  step  n 
cn(n),dn(n)  —  the  recursive  variables,  at  time  step  n 
rni,tni,cni,dni  —  the  above  variables  for  a  single  time 
step 

xin,  yin  —  the  action  and  angle  variables  at  the  current 

time  step,  passed  to  subroutine  map 
xout,  yout  —  the  variables  at  the  next  time  step,  output 

from  subroutine  map 
abnd  —  The  bound  on  en  and  dn,  from  the  theorem 
en(n)  —  the  same  as  dn(n) ,  but  backwards 
delta  —  the  noise  level 
E  —  the  energy 

xl,yl,xdl,ydl  —  the  phase  space  variable  initial 
conditions 

dimen  —  the  dimension  of  the  system 

uone(n) ,utwo(n)  —  two  choices  for  the  unstable  manifold, 

which  must  converge 
sone(n) ,stwo(n)  —  two  choices  for  the  stable  manifold 
udiff l,udiff2  —  the  convergence  criterion  for  the  manifolds 
ee  —  used  in  determining  of  e(n)  approaches  bound  of  theorem 
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*  1  —  flag  used  to  determine  end  of  shadowability 
* 

Double  Precision  x(S00001) ,y(500001) ,z,ajac(4,4) ,ajacinv(4,4) 

Double  Precision  sone(4, 500001), csc(lOOOO) ,csec 

Double  Precision  ul(4),u2(4) ,sl(4) ,s2(4),rn(10000),tn(10000) 

Double  Precision  cn(10000) ,dn(10000),bnd,xin,yin,xout,yout 

Double  Precision  abnd, en (10000) ,cscl , csc2,rni,tni ,ci ,di 

Double  Precision  bbnd,a,b,xl,yl 

Double  Precision  uone(4, 500000) 

Double  Precision  xdin,xdout,ydin,ydout,xdot( 500000) ,ydot(500000) 

Double  Precision  xdl ,ydl ,x2,y2,xd2,yd2 

Double  Precision  t ,xi,yi,xdi,ydi,h 

Double  Precision  xla(4) ,dxdt(4) ,dmdt(4) 

Double  Precision  ulnorm,u2norm,slnorm,s2norm 

Double  Precision  anorm(4) ,aident(4,4) 

Double  Precision  theta.sigma 

double  precision  delta, E.dimen 

Real  ran2 

Integer  jran 

Integer  I,i,k,l,kk,ltime,lreg 

data  idum/-3303051/ 

common  //  delta, E.dimen 

1=5000 

lreg=0 

ltime=0 

bnd=2.d0 

cn(l)=0.d0 

abnd=dimen**(-5.d0/2.d0)*bnd**(-2.d0)/dSqrt (delta) 

dn(I)=abnd 

en(l)=dn(I) 

1=0 

ii=0 

h=l.d-4 

do  i=l,4 

do  J»l,4 
aident(i,j)=0.d0 
end  do 
end  do 
do  i=l,4 

aident(i,i)=l.d0 
end  do 

read(21,*)i,sl(l),sl(2),sl(3),sl(4),lreg 
read(22,^)i,ul(l) ,ul(2) ,ul(3) ,ul(4) ,lreg 
read(23,*)xl,yl,xdl,ydl 
if  ((lreg  .eq.   1)  .or.   (lreg  .eq.2))  then 

return 
endif 

write (6,*) 'stable  mild' 
write(6,*)xl,yl,xdl,ydl 
Hrite(6,*)sl(l).sl(2),sl(3),sl(4) 
do  j=l,4 

uone(j,l)=ul(j) 

sone(j ,5001)=sl(j) 
end  do 
x(l)=xl 
y(D=yi 
xdot(l)=xdl 
ydot(l)=ydl 
* 

*  Begin  sain  loop  to  calculate  shadowing  times,  locations. 

Do  kk= 1,20000 

*  Calculate  the  unstable  manifold  for  the  next  I  time  steps. 

Do  i=l,I 

xi=x(i) 
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yi=y(i) 
xdi=xdot(i) 
ydi=ydot(i) 
anorm (  1 )  =uone  ( 1 ,  i) 
anorm ( 2 ) =uone ( 2 , i ) 
anorm ( 3 ) =uone ( 3 , i ) 
anorm ( 4 ) =uone (4 , i) 
Call  Jacob(ajac,zi,yi) 
do  k=l,4 
do  j-1,4 

a j  ac (k , j ) =a j  ac (k , j ) *h+aident (k , j ) 
end  do 

anorm(k)=ajac(k,l)*ul(l)+ajac(k,2)*ul(2)+ 
!  ajac(k,3)*ul(3)+ajac(k,4)*ul(4) 

end  do 

Call  Normal (anorm) 
uone(l,i+l)=anorm(l) 
uone(2,i+l)=anorm(2) 
uone(3,i+l)=anorm(3) 
uone(4,i+l)=anorm(4) 
anorm(l)=x(i) 
anorm(2)=y(i) 
anorm ( 3 ) =xdot ( i ) 
anorm(4)=ydot(i) 
call  Taylor ( anorm, h) 
x(i+l)=anorm(l) 
y(i+l)=anorm(2) 
xdot(i+l)=anorm(3) 
ydot(i+l)=anorm(4) 

if  ((dabs(x(i+l))  .gt.   I.d03)  .or.   (dabs(y(i+l))  .gt.   I.d03)) 
then 

lreg  =  2 
return 
endif 
End  do 
* 

*  Calculate  the  stable  manifold  for  the  next  K  time  steps. 
* 

Do  i=I,l,-l 
xi=x(i) 
yi=y(i) 
xdi=xdot(i) 
ydi=ydot(i) 
anorm(l)=sone(l,i+l) 
anorm (2 ) =sone ( 2 , i+1) 
anorm(3)=sone(3,i+l) 
anorm (4 ) = s  one ( 4 , i+1 ) 
Call  Jacinv(ajacinv,xi,yi) 
do  k=l,4 
do  j=1.4 

ajacinv(k, j)=ajacinv(k, j)*h+aident(k, j) 
end  do 

anorm(k)=ajacinv(k,l)*sl(l)+ajacinv(k,2)*sl(2)+ 
!  ajacinv(k,3)*sl(3)+ajacinv(k,4)*sl(4) 

end  do 

Call  lormal(anorm) 
sone ( 1 , i) = anorm (1 ) 
sone(2,i)=anorm(2) 
sone(3,i)=anorm(3) 
sone (4 , i) =anorm(4) 
End  do 
* 

*  Given  the  stable  and  unstable  manifolds,  determine  the  important 

*  shadowing  quantities. 


* 


Do  i=l,I 
do  j=l,4 
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ul(j)=uone(j,i) 

sl(j)=sone(j,i) 

and  do 
Call  angle(ul ,sl , csec.theta) 
csc(i)=csec 
do  j-i.4 

ul(j)=uone(j ,i+l) 

sl(j)=sone(j,i+l) 

end  do 
Call  angle  (ul.sl ,csec,theta) 
csc(i+lj=csec 
xin=x(i) 
yin=y(i) 
do  j-1.4 

ul(j)=uone(j,i+l) 

»l(j)=sone(j,i) 

end  do 
xi=x(i-l) 
yi=y(i-l) 
xdi=xdot(i-l) 
ydi=ydot(i-l) 
Call  Jacob(ajac,xi,yi) 
Call  Jacinv(ajacinv,xi,yi) 

Call  bound(ajac,ajacinv,ul,sl,rni,tni,h,aident) 
rn(i)=rni 
tn(i)=tni 
ci=cn(i) 
di=en(i) 
C8cl=csc(i) 
C8c2=csc(i+1) 

Call  recur (cscl , cs c2 , rni , tni , ci , di , abnd , 1 , It ime ) 
if  (1  .eq.   1)  then 

x2=x(i-l) 

y2=y(i-l) 

xd2=xdot(i-l) 

yd2=ydot(i-l) 

do  J-1.4 

anor»(j)=uone(j ,i-l) 

end  do 

return 
end  if 
cn(i+l)=ci 
en(i+l)=di 
End  do 
* 

*  If  glitch  has  not  yet  been  encountered,  re-do  process  for  the  next 

*  I  time  steps. 

x(l)=x(I+l) 

y(l)*y(I+l) 

xdot(l)=xdot(I+l) 

ydot(l)=ydot(I+l) 

en(l)=en(I+l) 

cn(l)=cn(I+l) 

anor*(l)=uone(l,I+l) 

anor«(2)=uon«(2,I+l) 

anorm(3)  =uone  (3 ,1+1 ) 

anorm ( 4 ) =uone ( 4 , N+ 1 ) 

Call  normal (anorm) 

uone (1,1) =  anorm ( 1 ) 

uone(2,l)=anora(2) 

uone (3 , 1 ) = anorm ( 3 ) 

uone(4,l)=anora(4) 

xl=x(l) 

yi=y(D 

xdl=xdot(l) 

ydl=ydot(l) 
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ii=20000 

call  smfld(xl,yl,xdl,ydl,sl,lreg,ii) 

if  ((lreg  . eq.   1)  .or.   (lreg  . eq.2))  then 

return 
endif 

8one(l,I+l)=sl(l) 
sone(2,I+l)=8l(2) 
sone(3,I+l)=sl(3) 
sone(4.I+l)=sl(4) 
end  do 

Write(6,*)  'map  can  be  shadowed  lor  over  100  Billion  iterates.' 
End 

Subroutine  Normal (anorm) 

*  The  four-dimensional  vector  anorm  is  normalized. 

Double  Precision  anorm(4),w 

w=DSqrt(anorm(l)**2+anorm(2)**2+anorm(3)**2+anorm(4)**2) 

anorm(l)=anorm(l)/v 

anorm(2)=anorm(2)/w 

anorm(3)=anorm(3)/w 

anorm(4)=anorm(4)/B 

Return 

End 

Subroutine  angle(u, s ,csec,theta) 

*  The  angle  between  the  stable  and  unstable  manifolds  at  a  given  time 

*  step  is  calucated  for  one-dimensional  manifolds  in  a  four-dimensional 

*  space . 

Double  Precision  u(4) ,s(4) ,csec 

Double  Precision  theta, dot ,normu, norms, alpha 

dot=0.d0 

normu-0 . dO 

norms =0.d0 

do  i=l,4 

dot=dot+u(i)*s(i) 

normu=normu+u ( i ) *  *2 

norms =norms + s ( i ) * * 2 
end  do 

normu=Dsqrt (normu) 
norms=Dsqrt (norms) 
alpha=dot / ( no  r mu  *  no  r ms ) 
if  (alpha  .gt.   l.dO)  then  alpha  =  l.dO 
theta=dacos (alpha) 
if  (Ds in (theta;  .gt.   0.d0)then 
csec=Dabs(l.dO/Ds in (theta)) 
endif 
Return 
End 

Subroutine  bound (aj ac , aj acinv, u, s , rn , tn , delt .aident) 

*  The  bound  on  the  growth  of  one-dimensional  stable  and  unstable  manifolds 

*  is  calculated. 

Double  precision  ajac(4,4),ajacinv(4,4) ,u(4) ,s(4) 
Double  precision  sj(4) ,uj(4) , delt, aident (4, 4) 
Double  precision  rn,tn,sigmal,sigma2 
do  i=l,4 

sj(i)=0.d0 
uj(i)=0.d0 
end  do 
do  i=l,4 
do  j-1,4 

sj(i)=sj(i)+(ajac(i,j)*delt+aident(i,j))*s(j) 
uj(i)=uj(i)+(ajacinv(i,j)*delt+aident(i,j))*u(j) 
end  do 
end  do 
rn=Dsqrt(sj (l)**2+sj (2)**2+s j (3)**2+s j (4)**2) 
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tn=Dsqrt(uj (l)**2+uj (2)**2+uj (3)**2+uj (4)**2) 
return 

end 

Subroutine  recur (cscl , csc2 ,rni ,tni , ci , ei ,abnd, 1 , ltime) 

*  This  subroutine  calculates  the  resursive  variable  ci  and  ei  for  a 

*  four-dimensional  system,  given  the  angle  between  the  stable  and  unstable 

*  manifolds,  the  bounds  on  the  growth  of  these  manifolds,  and  the 

*  bound  of  the  shadowing  theorem.    The  time  step  is  increased  during  this 

*  subroutine.    If  the  shadowing  theorem  bound  is  surpassed,  a  flag  is 

*  sent  to  the  main  program. 

Double  precision  cscl,csc2,abnd,en 

Double  precision  mi,tni,ci,ei 

ltime=ltime+l 

ci=csc2+rni*ci 

en=(ei-cscl)/tni 

if  (en  .gt.   abnd)  then 

ei=abnd 

else 

ei=en 

endif 

if  ((ci  .gt.   abnd)  .or.   (ei  .It.  0.0))  then 

write (6,*)  'unstable  point  is  no  longer  shadowable' 

write (6,*)  ltime 

vrite(6,*)'cscl,csc2' ,cscl,csc2 

write(6 , *) 'rn.tn' ,rni,tni 

write(6,*)'ci,ei' ,ci,ei 

1=1 

endif 

return 

End 

Subroutine  Jacob(ajac,xin,yin) 

Double  Precision  ajac(4,4),xin,yin 

if  (xin  .It.   l.d-20)  then 
xin=0 . dO 

endif 

if  (yin  .It.   l.d-20)  then 
yin=0.d0 

endif 

ajac(l,l)=0.d0 

ajac(l,2)=0.d0 

ajac(l,3)=l.d0 

ajac(l,4)=0.d0 

ajac(2,l)=0.d0 

ajac(2,2)=0.d0 

ajac(2,3)=0.d0 

ajac(2,4)=l.d0 

ajac(3,l)=  -l.dO  -  6.d0*xin**2  -  6.d0*xin**4  -  2.d0*yin  - 
!  12.d0*xin**2*yin  -  2.d0*yin**2  - 

!  12.d0*xin**2*yin**2  -  4.d0*yin**3/3.d0  - 

!  2.d0*yin**4/3.d0 

ajac(3,2)=-2.d0*xin  -  4.d0*xin**3  -  4.d0*xin*yin  - 
!  8.d0*xin**3*yin  -   4.d0*xin*yin**2  - 

!  8.d0*xin*yin**3/3.d0 

ajac(3,3)=0.d0 

ajac(3,4)=0.d0 

ajac(4,l)=-2.d0*xin  -  4.d0*xin**3  -  4.d0*xin*yin  - 
!  8.d0*xin**3*yin  -   4.d0*xin*yin**2  - 

!  8.d0*xin*yin**3/3.d0 

ajac(4,2)=-l.d0  -  2.d0*xin**2  -  2.d0*xin**4  + 

2.d0*yin  -  4.d0*xin**2*yin  -  6.d0*yin**2  - 
4.d0*xin**2*yin**2  +  20.d0*yin**3/3.d0  - 
22.d0*yin**4/3.d0 

ajac(4,3)=0.d0 

ajac(4,4)=0.d0 

return 
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end 


Subroutine  Jacinv(ajacinv,xin,yin) 

Double  Precision  ajacinv(4,4) ,xin,yin,denom 

if  (xin  .It.   l.d-20)  then 
x  in=0 . dO 

endiz 

if  (yin  .It.   l.d-20)  then 
yin=0 . dO 

endif 

denoB=(l.dO  +  4.d0*xin**2  +  4.d0*xin**4  +  8.d0*xin**6  + 
!        12.d0*xin**8  -  8.d0*xin**2*yin  -  24.d0*xin**4*yin  - 
!        16 . d0*xin**6*yin  +  4.d0*yin**2  +  8.d0*xin**2*yin**2  + 
!        40.d0*xin**4*yin**2  -  16.d0*xin**6*yin**2  +  8.d0*yin**3/3.d0 
!        -  16.d0*xin**2*yin**3  -  80.d0*xin**4*yin**3/3.d0  + 
!        4.d0*yin**4  +  40.d0*xin**2*yin**4/3.d0  + 
!        152.d0*xin**4*yin**4/3.d0  +  8.d0*yin**5  - 
!        16.d0*xin**2*yin**5/3.d0  +  88.d0*yin**6/9.d0  + 
!        752.d0*xin**2*yin**6/9.d0  +  16.d0*yin**7/3.d0  + 
!        44.d0*yin**8/9.d0) 

ajacinv(l,l)=0.dO 

ajacinv(l,2)=0.d0 

ajacinv(l,3)=(-l.d0  -  2.d0*xin**2  -  2.d0*xin**4  +  2.d0*yin  - 
!  4.d0*xin**2*yin  -  6.d0*yin**2  -  4.d0*xin**2*yin**2  + 
!        20.d0*yin**3/3.d0  -  22.d0*yin**4/3.d0)/denoa 

ajacinv(l,4)=(2.d0*xin  ♦  4.d0*xin**3  +  4.d0*xin*yin  + 
!  "      8.d0*xin**3*yin  +  4.d0*xin*yin**2  + 
!        8.d0*xin*yin**3/3.d0)/denom 

ajacinv(2,l)=0.d0 

ajacinv(2,2)-0.d0 

ajacinv(2,3)=ajacinv(l ,4) 

ajacinv(2,4)=(-l.d0  -  6.d0*xin**2  -  6.d0*xin**4  -  2.d0*yin  - 
!  "      12.d0*xin**2*yin  -  2.d0*yin**2  -  12.d0*xin**2*yin**2  - 
!        4.d0*yin**3/3.d0  -  2.d0*yin**4/3.d0)/denom 

ajacinv(3,l)=l.d0 

ajacinv(3,2)=0.d0 

ajacinv(3,3)=0.d0 

ajacinv(3,4)=0.d0 

ajacinv(4,l)=0.d0 

ajacinv(4,2)=l.d0 

ajacinv(4,3)=0.d0 

ajacinv(4,4)=0.d0 

return 

end 

Subroutine  Taylor(am,h) 

*  This  subroutine  uses  the  Taylor  method  to 

*  solve  the  Sixth  order  truncation  of  the  Toda  lattice  potential. 

double  precision  am(4) ,x(0:10) ,y(0:10),xh,xhd,yh,yhd,h,t 
double  precision  xl,yl,x2,y2 

*  Set  the  initial  conditions 

k=5 

x(0)=am(l) 
y(0)=a«(2) 
x(l)=aM(3) 
y(l)=a«(4) 
* 

x(2)=-x(0)  -  2.d0*x(0)**3  -  6.d0*x(0)**5/5.d0  -  2.d0*x(0)*y(0)  - 

4.d0*x(0)**3*y(0)  -  2.d0*x(0)*y(0)**2  - 

4.d0*x(0)**3*y(0)**2  -  4.d0*x(0)*y(0)**3/3.d0  - 

2.d0*x(0)*y(0)**4/3.d0 
y(2)=-x(0)**2  -  x(0)**4  -  y(0)  -  2.d0*x(0)**2*y(0)  - 

2.d0*x(0)**4*y(0)  +  y(0)**2  -  2.d0*x(0)**2*y(0)**2  - 

2.d0*y(0)**3  -  4.d0*x(0)**2*y(0)**3/3.d0  +  5.d0*y(0)**4/3.d0  - 

22.d0*y(0)**5/15.d0 
x(3)=-x(l)  -  6.d0*x(0)**2*x(l)  -  6.d0*x(0)**4*x(l)  - 

2.d0*y(0)*x(l)  - 
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12.d0*x(0 
12.d0*x(0 
2.d0*y(0) 

2.d0*x(0) 
8.d0*x(0) 

8.d0*x(0) 


y(3)=-2.d0*x(0)*x(l)  -  4.d0*x(0)**3*x(l)  -  4.dO*x(0)*y(0)*x(l) 


8.d0*x(0) 

8.d0*x(0) 

2.d0*x(0) 

6.d0*y(0) 

2O.dO*y(0 
x(4)=-12.d0*x 

24.d0*x(0 

4.d0*x(l) 

24.d0*x(0 

48.d0*x(0 

8.d0*y(0) 

4.dO*x(0) 

8.d0*x(0) 

8.d0*x(0) 

6.d0*x(0) 

2.d0*y(0) 

4.dO*y(0) 

2.d0*y(0) 

4.dO*x(0) 

4.dO*x(0) 
y(4)=-2.d0*x( 

24.d0*x(0 

8.d0*y(0) 

16.d0*x(0 

16.d0*x(0 

4.d0*x(0) 

8.dO*x(0) 

20.d0*y(0 

2.dO*x(0) 

4.d0*x(0) 

8.d0*x(0) 

4.d0*x(0) 

2.d0*x(0) 

4.d0*x(0) 

4.d0*x(0) 

20.d0*y(0 
xl=-12.d0*x(l 

24.d0*y(0 

72.d0*x(0 

12.d0*x(l 

24.d0*y(0 

8.dO*x(0) 

36.d0*x(0 

72.d0*x(0 

6.d0*y(l) 

12.d0*y(0 

12.d0*y(0 
x2=-36.d0*x(0 

72.d0*x(0 

8.d0*y(0) 

24.d0*x(0 

24.d0*x(0 

6.d0*x(0) 

2.d0*y(0) 

4.d0*y(0) 

2.d0*x(0) 

8.d0*x(0) 

8.d0*x(0) 
x(5)=xl+x2 
yl=-24.d0*x(0)*x(l)**3 


**2*y(0)*x(l)  -  2.d0*y(0)**2*x(l)  - 

**2*y(0)**2*x(l)  -  4.d0*y(0)**3*x(l)/3.d0  - 

*4*x(l)/3.d0  - 

y(l)  -  4.d0*x(0)**3*y(l)  -  4.d0*x(0)*y(0)*y(l) 

*3*y(0)*y(l)  -  4.d0*x(0)*y(0)**2*y(l)  - 

y(0)**3*y(l)/3.d0 


*3*y(0)*x(l)  -  4.d0*x(0)*y(0)**2*x(l) 

y(0)**3*x(l)/3.d0  -  y(l)  -  2.d0*x(0)**2*y(l)  - 

*4*y(l)  +  2.d0*y(0)*y(l)  -  4.d0*x(0)**2*y(0)*y(l)  - 

*2*y(l)  -  4.d0*x(0)**2*y(0)**2*y(l)  + 

**3*y(l)/3.d0  -  22.d0*y(0)**4*y(l)/3.d0 

0)*x(l)**2  -  24.d0*x(0)**3*x(l)**2  - 

*y(0)*x(l)**2  -  24.d0*x(0)*y(0)**2*x(l)**2  - 

y(l)  - 

**2*x(l)*y(l)  -  8.d0*y(0)*x(l)*y(l)  - 

♦*2«y(0)*x(l)*y(l)  - 

*2*x(l)*y(l)  -  16.d0*y(0)**3*x(l)*y(l)/3.d0  - 

y(l)**2  -  8.d0*x(0)**3*y(l)**2  - 

y(0)*y(l)**2  - 

y(0)**2*y(l)**2  -  x(2)  -  6.d0*x(0)**2*x(2)  - 

*4*x(2)  -  2.d0*y(0)*x(2)  -  12.d0*x(0)**2*y(0)*x(2) 

*2*x(2)  -   12.d0*x(0)**2*y(0)**2*x(2)  - 

*3*x(2)/3.d0  - 

*4*x(2)/3.d0  -  2.dO*x(0)*y(2)  -  4.d0*x(0)**3*y(2)  - 

y(0)*y(2)  -  8.d0*x(0)**3*y(0)*y(2)  - 

y(0)**2*y(2)  -  8.d0*x(0)*y(0)**3*y(2)/3.d0 

)**2  -  12.d0*x(0)**2*x(l)**2  -  4.d0*y(0)*x(l)**2  - 

**2*y(0)*x(l)**2  -  4.d0*y(0)**2*x(l)**2  - 

*3*x(l)**2/3.d0  -  8.d0*x(0)*x(l)*y(l)  - 

**3*x(l)*y(l)  -  16.d0*x(0)*y(0)*x(l)*y(l)  - 

*y(0)**2*x(l)*y(l)  +  2.d0*y(l)**2  - 

♦2*y(l)**2  -  12.d0*y(0)*y(l)**2  - 

*2*y(0)*y(l)**2  + 

**2*y(l)**2  -  88.d0*y(0)**3*y(l)**2/3.d0  - 

x(2)  - 

*3«x(2)  -  4.d0*x(0)*y(0)*x(2)  - 

*3*y(0)*x(2)  - 

y(0)**2*x(2)  -     8.d0*x(0)*y(0)**3*x(2)/3.d0  -  y(2) 

*2*y(2)  -  2.d0*x(0)**4*y(2)  +  2.d0*y(0)*y(2)  - 

*2*y(0)*y(2)  -  6.d0*y(0)**2*y(2)  - 

*2*y(0)**2*y(2)  + 

**3*y(2)/3.d0  -  22.d0*y(0)**4*y(2)/3.d0 

**3  -  72.d0*x(0)**2*x(l)**3  - 

*x(l)**3  -  24.d0*y(0)**2*x(l)**3  - 

*x(l)**2*y(l)  -   144.d0*x(0)*y(0)*x(l)**2*y(l)  - 

*y(l)**2  -  72.d0*x(0)**2*x(l)*y(l)**2  - 

*x(l)*y(l)**2  -  24.dO*y(0)**2*x(l)*y(l)**2  - 

y(l)**3  -   16.d0*x(0)*y(0)*y(l)**3  - 

*x(l)*x(2)  -  72.d0*x(0)**3*x(l)*x(2)  - 

*y(0)*x(l)*x(2)  -  72.d0*x(0)*y(0)**2*x(l)*x(2)  - 

x(2)  -  36.d0*x(0)**2*y(l)*x(2)  -  6.d0*x(l)*y(2)  - 

*y(l)*x(2)  -  72.d0*x(0)**2*y(0)*y(l)*x(2)  - 

**2*y(l)*x(2)  -  8.d0*y(0)**3*y(l)*x(2) 

**2*x(l)*y(2)  -   12.d0*y(0)*x(l)*y(2)  - 

**2*y(0)*x(l)*y(2)  -   12.d0*y(0)**2*x(l)*y(2)  - 

*3*x(l)*y(2)  -   12.d0*x(0)*y(l)*y(2)  - 

**3*y(l)*y(2)  -  24.d0*x(0)*y(0)*y(l)*y(2)  - 

*y(0)**2*y(l)*y(2)  -  x(3)  -  6.d0*x(0)**2*x(3)  - 

♦4*x(3)  -  2.d0*y(0)*x(3)  -   12.d0*x(0)**2*y(0)*x(3) 

*2*x(3)  -   12.d0*x(0)**2*y(0)**2*x(3)  - 

*3*x(3)/3.d0  -  2.d0*y(0)**4*x(3)/3.d0  - 

y(3)  -  4.d0*x(0)**3*y(3)  -  4.d0*x(0)*y(0)*y(3)  - 

*3*y(0)*y(3)  -  4.d0*x(0)*y(0)**2*y(3)  - 

y(0)**3*y(3)/3.d0 

48.d0*x(0)*y(0)*x(l)**3  - 
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12.d0*x(i)**2*y(l)  - 

72.d0*x(0)**2*x(l)**2*y(l)  -  24.d0*y(0)*x(l)**2*y(l)  - 

24.d0*y(0)**2*x(l)**2*y(l)  -  24.d0*x(0)*x(l)*y(l)**2  - 

48.d0*x(0)*y(0)*x(l)*y(l)**2  -  12.d0*y(l)**3  - 

8.d0*x(0)**2*y(l)**3  + 

40.d0*y(0)*y(l)**3  -  88.d0*y(0)**2*y(l)**3  -  6.d0*x(l)*x(2)  - 

36.d0*x(0)**2*x(l)*x(2)  -  12.d0*y(0)*x(l)*x(2)  - 

72.d0*x(0)**2*y(0)*x(l)*x(2)  -  12.d0*y(0)**2*x(l)*x(2)  - 

8.d0*y(0)**3*x(l)*x(2)  -  12.d0*x(0)*y(l)*x(2)  - 

24.d0*x(0)**3*y(l)*x(2)  - 

24.d0*x(0)*y(0)*y(l)*x(2)  -  24.d0*x(0)*y(0)**2*y(l)*x(2)  - 

12.dO*x(0)*x(l)*y(2)  -  24.d0*x(0)**3*x(l)*y(2) 

y2=-24.d0*x(0)*y(0)*x(l)*y(2)  - 

24.d0*x(0)*y(0)**2*x(l)*y(2)  +  6.d0*y(l)*y(2)  - 
12.d0*x(0)**2*y(l)*y(2)  -  36.d0*y(0)*y(l)*y(2)  - 
24.d0*x(0)**2*y(0)*y(l)*y(2)  +  60.d0*y(0)**2*y(l)*y(2)  - 
88.d0*y(0)**3*y(l)*y(2)  -  2.d0*x(0)*x(3)  -  4.d0*x(0)**3*x(3) 
4.d0*x(0)*y(0)*x(3)  -  8.d0*x(0)**3*y(0)*x(3)  - 
4d0*x(0)*y(0)**2*x(3)  - 

8.d0*x(0)*y(0)**3*x(3)/3.d0  -  y(3)  -  2.d0*x(0)**2*y(3)  - 
2.d0*x(0)**4*y(3)  +  2.d0*y(0)*y(3)  -  4.d0*x(0)**2*y(0)*y(3)  - 
6.d0*y(0)**2*y(3)  -  4.d0*x(0)**2*y(0)**2*y(3)  + 
20.d0*y(0)**3*y(3)/3.d0  -  22.d0*y(0)**4*y(3)/3.d0 

y(5)=yl+y2 

xh=0.dO 
xhd=0.dO 
yh-O.dO 
yhd=0.dO 
do  i=0,k 

xh=xh+x ( i) *h**i/iactrl ( i) 

xhd=xhd+x(i+l)*h**i/iactrl(i) 

yh=yh+y(i)*h**i/iactrl(i) 

yhd=yhd+y(i+l)*h**i/iactrl(i) 

end  do 
aa(l)=xh 
am(2)=yh 
am(3)=xhd 
an(4)=yhd 
•nd 

Subroutine  derivs(t ,x,dxdt) 

double  precision  x(4) ,dxdt(4) ,t 

dxdt(l)=x(3) 

dxdt(2)=x(4) 

dxdt(3)=-x(l)  -  2.d0*x(l)**3  -  6.dO*x(l)**S/5.dO  - 
!    2.d0*x(l)*x(2)  -  4.d0*x(l)**3*x(2)  -  2.d0*x(l)*x(2)**2  - 
!    4.d0*x(l)**3*x(2)**2  -  4.d0*x(l)*x(2)**3/3.d0  - 
!    2.d0*x(l)*x(2)**4/3.d0 

dxdt(4)=-x(l)**2  -  x(l)**4  -  x(2)  -  2.d0*x(l)**2*x(2)  - 
!    2.d0*x(l)**4*x(2)  +  x(2)**2  -  2.d0*x(l)**2*x(2)**2  - 
!    2.d0*x(2)**3  -  4.d0*x(l)**2*x(2)**3/3.d0  +  S.d0*x(2)**4/3.d0 
!    22.d0*x(2)**5/15.d0 

return 

end 


Subroutine  u«ild(xl,yl,xdl,ydl,ul,lreg,ii) 

Double  Precision  x(5002) ,y(5002) ,z,ajac(4,4) ,ajacinv(4,4) 

Double  Precision  ul(4),u2(4) 

Double  Precision  xl ,yl .xdl.ydl 

Double  Precision  xd(5002) ,yd(5002) 

Double  Precision  xi.yi.h 

Double  Precision  anorm(4) ,adent(4,4) 

Double  Precision  theta,csec,dmdt(4) ,t 

double  precision  delta, E.dimen 
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Integer  i.k.l.kk.lreg 
common  //  delta.E.dimen 
t=0.dO 
do  i=l,4 
do  j=1.4 

adent(i, j)=0.dO 
end  do 
end  do 
do  i=l,4 

adent(i,i)=l.dO 
end  do 
lreg=0 
ltiae=0 
ii=0 
iii=0 
x(l)=xl 
y(l)=yi 
zd(l)=zdl 
yd(l)=ydl 
* 

*  Set  initial  choices  for  the  unstable  manifold 
* 

anorm(l)=l.d0/3.d0 
anorm(2)=l.d0 
anorm(3)=l.d0 
anorm(4)=l.d0 
Call  Mormal (anorm) 
ul(l)=anorm(l) 
ul(2)=anorm(2) 
ul(3)=anorm(3) 
ul(4)=anor«(4) 
anorm(l)=l.dO 
anorm(2)=l.d0/3.d0 
anor*(3)=l.d0 
anorm(4)=l.d0 
Call  Normal (anorm) 
u2(l)=anorm(l) 
u2(2)=anorm(2) 
u2(3)=anorm(3) 
u2(4)=anorm(4) 
h=l.d-2 
* 

*  Make  sure  that  the  two  choices  for  the  unstable  manifold  converge. 

* 

i=l 
100  xi=x(i) 

yi=y(i) 

Call  Jacob(ajac,xi,yi) 
do  k=l,4 
do  j-1,4 

ajac(k,j)=ajac(k,j)*h+adent(k,j) 
end  do 

anorm(k)=ajac(k,l)*ul(l)+ajac(k,2)*ul(2)+ 
!  ajac(k,3)*ul(3)+ajac(k,4)*ul(4) 

end  do 

Call  Mormal (anorm) 
ul(l)=anorm(l) 
ul(2)=anorm(2) 
ul(3)=anorm(3) 
ul(4)=anorm(4) 
do  k=1.4 

anorm(k)=ajac(k,l)*u2(l)+ajac(k,2)*u2(2)+ 
!  ajac(k,3)*u2(3)+ajac(k,4)*u2(4) 

end  do 

Call  Normal (anorm) 

u2(l)=anorm(l) 

u2(2)=anorm(2) 
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u2(3)=anorm(3) 
u2(4)=anorm(4) 
anorm(l)=x(i) 
anora(2)=y(i) 
anora(3)=xd(i) 
anorm(4)=yd(i) 
call  Taylor (anorm.h) 
x(i+l)=anorm(l) 
y(i+l)=anorm(2) 
xd(i+l)=anorm(3) 
yd(i+l)=anorm(4) 

if  ((dabs(x(i+l))  .gt.   I.d03)  .or. 
!  (dabs(y(i+l))  .gt.   I.d03))  then 

lreg  =  2 
write(19,*)x(i+l) ,y(i+l) 
return 
endif 

call  angle (ul,u2,csec,theta) 
if  (Dsin(theta)  .It.   l.d-6)then 
ii=i 
xl=x(i) 
yl=y(i) 
xdl=xd(i) 
ydl=yd(i) 
return 
endif 
i=i+l 

iz  (i  .gt.   5000)  then 
lreg  =  1 
return 
endif 
gotolOO 
end 


Subroutine  smfld(xl,yl,xdl,ydl,sl,lreg,ii) 
Double  Precision  x(20001) ,y(20001) ,z,ajac(4,4) ,ajacinv(4,4) 
Double  Precision  sl(4),s2(4) 
Double  Precision  xl.yl ,xdl ,ydl 
Double  Precision  xd(20001) ,yd(20001) 
Double  Precision  xi.yi.h 
Double  Precision  anorm(4),adent(4,4) 
Double  Precision  theta,csec,dmdt(4) 
double  precision  delta, E.dimen 
Integer  i.k.l.kk.lreg 
conunon  //  delta, E.dimen 
do  i=l,4 
do  j=l,4 

adent(i,j)=0.d0 
end  do 
end  do 
do  i=l,4 

adent(i,i)=l.dO 
end  do 

if  ((lreg  .eq.   1)  .or.   (lreg  .eq.2))  then 
return 
endif 
lreg=0 
iii=0 
x(ii)=xl 
y(ii)=yl 
xd(ii)=xdl 
yd(ii)=ydl 
do  i=ii,ii+5000 

anor«(l)=x(i) 

anorm(2)=y(i) 

anorm(3)=xd(i) 
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anorm(4)=yd(i) 
h=l.d-4 

call  Taylor (anorm.h) 
x(i+l)=anorm(l) 
y(i+l)=anorm(2) 
xd(i+l)=anorm(3) 
yd(i+l)=anorm(4) 

if  ((dabs(x(i+l))  .gt.   I.d03)  .or. 
!  (dab»(y(i+l))  .gt.   I.d03))  then 

lreg  =  2 
write(19,*)x(i+l),y(i+l) 
return 
endif 
•nd  do 
do  i-ii-f 5000, 20000 
anorm(l)=x(i) 
anorm(2)=y(i) 
anorm(3)=xd(i) 
anorm(4)=yd(i) 
h=l.d-2 

call  Taylor (anorm.h) 
x(i+l)=anorm(l) 
y(i+l)=anorm(2) 
xd(i+l)=anorm(3) 
yd(i+l)=anorm(4) 

if  ((dabs(x(i+l))  .gt.   I.d03)  .or. 
!  (dabs(y(i+D)  .gt.   I.d03))  then 

lreg  =  2 
write(19,*)x(i+l)  ,y(i+l) 
return 
endif 
end  do 
* 

*  Approximate  the  stable  direction  by  some  random  vector. 
* 

i=20000 
250   anorm(l)=l.d0 

anorm(2)=l.d0 

anor»(3)=l.d0/3.d0 

anorm(4)=l .d0 

Call  Normal (anorm) 

sl(l)=anora(l) 

sl(2)=anorm(2) 

sl(3)=anorm(3) 

sl(4)=anorm(4) 

anorm(l)=l.d0 

anorm(2)=l.d0 

anorm(3)=l.d0 

anorm(4)=l.d0/3.d0 

Call  Normal (anorm) 

s2(l)=anorm(l) 

s2(2)=anorm(2) 

s2(3)=anorm(3) 

s2(4)=anorm(4) 
* 

*  Make  sure  that  the  two  choices  for  the  stable  manifold  converge. 

300         xi=x(i) 
yi=y(i) 
Call  Jacinv(ajacinv,xi,yi) 
do  k=l,4 
do  J-1,4 
ajacinv(k,j)=ajacinv(k,j)*h+adent(k,j) 
end  do 

anorm(k)=ajacinv(k,l)*sl(l)+ajacinv(k,2)*sl(2)+ 

I  ajacinv(k,3)*sl(3)+ajacinv(k,4)*sl(4) 

end  do 
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Call  Hormal (anorm) 
sl(l)sanorm(l) 
sl(2)=anorm(2) 
sl(3)=anorm(3) 
■l(4)=anorm(4) 
do  k=l,4 

anorm(k)=ajacinv(k,l)*s2(l)+ajacinv(k,2)*s2(2)+ 
!  ajacinv(k,3)*s2(3)+ajacinv(k,4)*s2(4) 

end  do 

Call  Honaal (anorm) 
s2(l)=anonn(l) 
•2(2)=anora(2) 
■2(3)=anorm(3) 
s2(4)=anora(4) 
call  angle (sl,s2,csec,theta) 
if  (Dsin(theta)  .It.   l.d-6)then 
iii=i 
goto  400 
endif 
i=i-l 

if  (i  . eq.   ii+5000)  then 
lreg=l 

write (6 , *) ' lreg= ' , lreg 
return 
endif 
goto  300 
400   do  i=iii, ii+5000, -1 
xi=x(i) 
yi=y(i) 
Call  Jacinv(ajacinv,xi,yi) 
do  k-1,4 
do  J-1.4 

ajacinv(k, j)=ajacinv(k,j)*h+adent(k, j) 
end  do 

anorm(k)=ajacinv(k, l)*sl(i)+ajacinv(k,2)*sl(2)+ 
!  ajacinv(k,3)*sl(3)+ajacinv(k,4)*sl(4) 

end  do 

Call  Mormal ( anorm) 
sl(l)=anora(l) 
»l(2)=anorB(2) 
»l(3)=anora(3) 
»l(4)=anor»(4) 
end  do 
xl=x(ii) 
yl=y(ii) 
xdl=xd(ii) 
ydl=yd(ii) 
Return 
end 


APPENDIX  B 
Fortran  Codes  -  four-dimensional  systems 


Program  Shadow 

Double  precision  xl.yl.randl 

Integer  i.ii.kk, jran.lr 

open  (unit=18,iile='short .dat ' ,status='new' ) 

open  (unit=19,iile='short2.dat ' ,status='new' ) 

jran=3037 

do  i=0,9 

Call  randl(jran) 

xl=randl(jran)*l.D-14 

Call  randl(jran) 

yl=randl(jran)*l.d-14 

write(18,*)xl,yl 

end  do 

rewind  (unit=18) 

do  ii=l,10 

kk=0 

1=0 

lr=l 

read(18,*)xl,yl 

call  ■«ld(xl.yl,xdl,ydl,ll.lr) 

srite(19,*)xl,yl,ll 

if  (lr  .eq.   1)  then 

write(16,*)  'Point  is  within  regular  region' 
endif 
if  (lr  .eq.   2)  then 

write(16,*)  'Particle  escapes' 
endif 
end  do 
end 

Subroutine  mild(xl,yl,xdl,ydl,ll,lr) 

Double  Precision  x(6000) ,y(6000) ,z, J(2,2), Jinv(2,2) ,m(2) 
Double  Precision  sone(2, 6000) ,stwo(2, 6000) ,csc(6000) ,csec 
Double  Precision  ul(2),u2(2) ,sl(2),s2(2) ,rn(6000) ,tn(6000) 
Double  Precision  cn(6000) ,dn(6000) ,bnd,xin,yin,xout,yout 
Double  Precision  abnd,en(6000) ,cscl,csc2,rndum,tndum,cdum,dduin 
Double  Precision  delta, bbnd, a, b,xl,yl,xdl,ydl,mdim,s(2, 6000) 
Double  Precision  uone(2, 6000) ,utwo(2, 6000) ,udiiil,udiff2 
Double  Precision  xdin,xdout,ydin,ydout,xdot(6000) ,ydot(6000) 
Integer  I.i.k.l.kk.ll.lr 
delta=1.0d-14 
a=1.4d0 
b=0.3d0 
1=6000 
mdim=2 . dO 
lr=0 
11=0 

bnd=0.d0 
x(l)=xl 
y(l)=yl 
xdot(l)=xdl 
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then 
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ydot(l)=ydl 
m(l)=dsqrt(2.d0) 
■(2)=l.d0 
Call  Normal (m) 
do  jj=1.2 

uone(jj,l)=«(jj) 
end  do 
«(l)=b/3.d0 
■(2)=l.d0 
Call  Hormal(m) 
do  jj=l,2 

utwo(jj,l)=«(jj) 
end  do 
cn(l)=0.dO 

abnd=mdi«**(-5.dO/2.dO)*4.dO**(-2.dO)/dSqrt(delta) 
dn(I)=abnd 
en(l)=dn(I) 
1=0 

Do  k= 1,20000 
Do  i=l,I 

xin=x(i) 

yin=y(i) 

Call  Jacob(J,a,b,xin,yin) 

do  jj-1,2 

■(jj)=J(jj,l)*uone(l,i)+J(jj,2)*uone(2,i) 
end  do 

Call  Normal (m) 

do  jj*l,2 

uone(jj,i+l)=*(jj) 

end  do 

do  jj-1,2 

■(jj)=J(jj,l)*utwo(l,i)+J(jj,2)*utwo(2,i) 

end  do 

do  jj-1,2 

■(jj)=«(jj)-uone(jj,i+l)*(uone(l,i+l)*«(l)+uone(2,i+l)*M(2)) 

end  do 
Call  Normal (m) 
do  jj=l,2 

utHo(jj,i+l)=»(jj) 
end  do 

xin=x(i) 

yin=y(i) 

xdin-xdot(i) 

ydin=ydot(i) 

call  mapCxin.yin.xdin.ydin.xout ,yout,xdout ,ydout ,a,b) 

x(i+l)=xout 

y(i+l)=yout 

xdot(i+l)=xdout 

ydot(i+l)=ydout 

if  ((dabs(x(i+l))  .gt.   I.d03)  .or.   (dabs(y(i+l))  .gt.   I.d03)) 

lr  =  2 
return 
endii 
bnd=dmaxl(bnd,Dab»(x(i)),Dabs(y(i))) 
End  do 

bbnd=2.d0*bnd 
■(l)=d«qrt(2.d0) 
■(2)=l.d0 
Call  Normal (m) 
do  jj=l,2 

8one(jj,I+l)=m(jj) 
end  do 
a(l)=l.d0 
«(2)=d»qrt(2.d0) 
Call  Normal (m) 
do  jj=l,2 
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stwo(jj,I+l)=m(jj) 
end  do 
Do  i-1,1,-1 

xin=x(i) 

yin=y(i) 

Call  Jacinv(Jinv,a,b,xin,yin) 

do  jj=l,2 

»(jj)=Jinv(jj,l)*sone(l,i+l)+Jinv(jj,2)*sone(2,i+l) 
end  do 
Call  Kormal (m) 

do  jj»i,a 

sone(jj,i)=«(jj) 
end  do 
do  ji=l,2 
■(jj)=Jinv(jj,l)*8two(l,i+l)+Jinv(jj,2)*8two(2,i+l) 

end  do 
do  14*1.3 

m(jj)=m(jj)-sone(jj,i)*(sone(l,i)*n(l)+sone(2,i)*«(2)) 
end  do 
Call  lormal(m) 
do  jj=l,2 

8two(jj,i)=»(jj) 
end  do 
End  do 
Do  i= 1,5000 
do  jj-1,2 

ul(jj)=uone(jj,i) 
u2(jj)=utwo(jj,i) 
8l(jj)=sone(jj,i) 

82(jj)=8t¥0(jj,i) 

end  do 

Call  angle(ul ,u2,sl ,s2,csec) 

C8c(i)=csec 

do  jj»1.2 

ul(jj)=uone(jj,i+l) 

sl(jj)=8one(jj,i+l) 

u2(jj)=utwo(jj.i+l) 

«2(jj)=8tHo(jj,i+l) 

end  do 
Call  angle  (ul ,u2,sl ,s2,csec) 
C8c(i+l;=c8ec 
xin=x(i) 
yin=y(i) 

Call  Jacob(J,a,b,xin,yin) 
Call  Jacinv(Jinv,a,b,xin,yin) 
do  JJ-1,2 

ul(jj)=uone(jj,i) 

u2(jj)=utwo(jj,i) 

sl(jj)=8one(jj,i) 

s2(jj)=stwo(jj,i) 

end  do 
Call  bound(J, Jinv.ul ,u2 >sl,s2,rndum,tndum) 
rn(i)=rndujn 
tn(i)=tndua 
cdua=cn(i) 
ddu«=en(i) 
cscl=c8c(i) 
C8c2=csc(i+1) 

Call  recur(cscl,csc2,rndum,tndum,cdum,ddum1abnd,k,kk,l,ll) 
if  (1  .eq.   1)  return 
cn(i+l)=cdun 
en(i+l)=ddun 
End  do 
x(l)=x(5001) 
y(l)=y(5001) 
xdot(l)=xdot(5001) 
ydot(l)=ydot(5001) 
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en(l)=en(5001) 

cn(l)=cn(5001) 

end  do 

Write(6,*)  'map  can  ba  shadowed  for  over  100  Billion  iterates.' 

End 

Subroutine  Hormal(m) 

Double  Precision  m(2),w 

w=DSqrt(«(l)**2+«(2)**2) 

■<l)=«(l)/s 

■(2)=»(2)/u 

Return 

End 

Subroutine  angle (au , au2 , as , as2 , csec) 

Double  Precision  au(2) ,as(2) , csec.atm.atn 

Double  Precision  an(2,2),bn(2,2) ,bnu(2,2),bns(2,2) 

Double  Precision  sn(2,2),un(2,2),id(2,2) 

double  Precision  csecl ,csec2,csecdii 

Double  Precision  a(2,2) ,w(2) ,ivl(2) ,fv2(2) 

n=2 

np=2 

do  ji=l,2 

an(jj,l)=as(jj) 

an(jj,2)=au(jj) 

end  do 
Call  ■atinv(an,bn,n,np,indx) 
do  jj=l,2 

do  jjj=1.2 

bns(jj,jjj)=0.dO 
bnu(jj,jjj)=0.dO 
end  do 

end  do 
do  jj=l,2 

bns(l,jj)=bn(l,jj) 

bnu(2,jj)=bn(2,jj) 

end  do 
do  ii=l,2 

do  jj=1.2 

sn(ii,jj)=as(ii)*bns(l,jj) 
un(ii,jj)=au(ii)*bnu(2,jj) 
end  do 

end  do 
do  ii=l,2 

do  jj=l,2 

id(ii,jj)=an(ii,jj)+un(ii,jj) 
end  do 
end  do 
do  ii=l,2 

id(ii,ii)=id(ii,ii)-l.dO 
end  do 
do  ii=l,2 
do  jj-1.2 

iz  (id(ii.ji)  .gt.   l.d-5)then 

write(6,*) 'error  caused  in  matrix  inversion' 
stop 
end  if 
end  do 
end  do 
do  ii=l,2 
do  ij=l,2 

a(ii,jj)=sn(l,ii)*sn(l,jj)+sn(2,ii)*sn(2,jj) 
end  do 
end  do 

call  tredl(np>n,a,v,fvl>fv2) 
call  tqlrat(n,w,fv2,ierr) 
if  (ierr  .ne.  0)  then 

write(6,*) 'error  caused  in  eigenvalue  routine' 
stop 
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end  if 

csecl=Dsqrt(max(w(l) ,w (2))) 
do  11*1,3 

do  ji=l,2 

a(ii,jj)=un(l,ii)*un(l,jj)+un(2,ii)*un(2,jj) 

end  do 
end  do 

call  tredl(np,n,a,w,ivl,fv2) 
call  tqlrat(n,w,fv2)ierr) 
if  (ierr  . ne.   0)  then 

write(6,*) 'error  caused  in  eigenvalue  routine' 

■top 
end  if 

csec2=Dsqrt(max(v(l),w(2))) 
csecdif =dabs(csecl-csec2) 
if  (csecdif  . gt.   l.d-6)then 

write(6,*) 'error  caused  in  matrix  inversion' 

stop 

endif 
csec=csecl 
Return 
End 

Subroutine  bound ( J , Jinv ,ul ,u2 , si , s2 , rndum, tndum) 
Double  precision  J(2,2) , Jinv(2,2) ,ul(2) ,u2(2) ,sl(2) ,s2(2) 
Double  precision  sj(2,2),uj(2,2) ,sjt(2,2),ujt(2,2) 
Double  precision  rndum, tndum, signal ,sigma2,s(2, 2) ,u(2,2) 
Double  precision  id(2,2) ,a(2,2) ,b(2,2) ,w(2) ,fvl(2) ,fv2(2) 
np=2 
n=2 

u2(l)=0.d0 
u2(2)=0.d0 
s2(l)=0.d0 
s2(2)=0.d0 
do  ii=l,2 

s(ii,l)=sl(ii) 

s(ii,2)=s2(ii) 

u(ii,l)=ul(ii) 

u(ii,2)=u2(ii) 

end  do 
do  ii=l,2 

do  jj-1,2 

sj(ii,jj)=J(ii,l)*s(l,jj)+J(ii,2)*s(2,ij) 
uj(ii, j j)=Jixvv(ii,l)*u(l, j j)+Jinv(ii,2;*u(2, jj) 
end  do 

end  do 
do  ii=l,2 

do  Jj-1,2 

sjt(ii,jj)=sj(l,ii)*sj(l,jj)+sj(2,ii)*sj(2,jj) 
ujt(ii,jj)=uj(l,ii)*uj(l,jj)+uj(2,ii)*uj(2,jj) 
end  do 

end  do 
call  tredl(np,n>sjt,v,fvl,fv2) 
call  tqlrat(nls,fv2,ierr) 
if  (ierr  .ne.  0)then 

write(6,*) 'error  caused  in  eigenvalues  in  subroutine  bound' 

stop 

end  if 
sigmal=Dsqrt(max(v(l),w(2))) 
do  ii=l,2 

do  jj=l,2 

id(ii,ij)=0.d0 
id(l,l)=l.d0 
end  do 

end  do 
do  ii=l,2 

do  jj=l,2 

sj(ii,jj)=s(l,ii)*s(l,jj)+s(2,ii)*s(2,jj)-id(ii,jj) 
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end  do 

end  do 
do  ii=l,2 

do  jj=l,2 

Sjt(ii,jj)=8j(l,ii)*8j(l,jj)+Bj(2,ii)*8j(2,jj) 
end  do 

end  do 
call  tredl(np,n,sjt,v,fvl,fv2) 
call  tqlrat(n,w,iv2,ierr) 
if  (ierr  .ne.   0)then 

write(6,*) 'error  caused  in  eigenvalues  in  subroutine  bound' 

stop 

end  if 
sigma2=Dsqrt(max(w(l) ,w(2))) 
rndum=sigmal/dsqrt ( 1 . d0-sigma2) 
call  tredl(np,n,ujt ,s,lvl,iv2) 
call  tqlrat(n,«,iv2,ierr) 
if  (ierr  .ne.   0)then 

write(6,*) 'error  caused  in  eigenvalues  in  subroutine  bound' 

stop 

end  if 
sigmal=Dsqrt(max(v(l),B(2))) 
do  ii=l,2 

do  jj-1,2 

uj(ii,jj)=u(l,ii)*u(l,jj)+u(2,ii)*u(2,jj)-id(ii,jj) 
end  do 

end  do 
do  ii=l,2 

do  jj=1.2 

ujt(ii,jj)=uj(l,ii)*uj(l,jj)+uj(2,ii)*uj(2,jj) 
end  do 

end  do 
call  tredl(np,n,ujt,v,fvl,fv2) 
call  tqlrat(n,w,iv2,ierr) 
if  (ierr  .ne.  0)then 

write(6,*) 'error  caused  in  eigenvalues  in  subroutine  bound' 

stop 

end  if 
sigma2=Dsqrt(max(w(l) ,w(2))) 
tndua=sigmal/dsqrt ( 1 . d0-sigma2) 
Return 
End 

Subroutine  recur ( cscl , csc2 , radum , tndum , cdum , edum , abnd , k , kk , 1 , 11 ) 
Double  precision  cscl, csc2, abnd 
Double  precision  radua, tndum, cdum, edum 
11=11+1 

cdum=csc2+rndum*cdum 
edum=dminl ( (edum-cscl) /tndum, abnd) 
if  ((cdum  .gt.  abnd)  .or.   (edum  .It.  0.0))  then 
write (6,*)  'unstable  point  is  no  longer  shadouable' 
write(6,*)  'length  is  approximately  ',k*5000,'  iterates.' 
write(6,*)  11 
1=1 

kk=k*5000 
endif 
return 
End 

Subroutine  Jacob(J,a,b,xin,yin) 
Double  Precision  J(2,2) ,a,b,xin,yin 
J(l,l)=-2.d0*xin 
J(1.2)=b 
J(2,l)=l.d0 
J(2,2)=0.d0 
return 
end 

Subroutine  Jacinv ( J inv , a ,b, xin ,yin) 
Double  Precision  Jinv(2,2),a,b,xin,yin 
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Jinv(l,l)=0.dO 

Jinv(l,2)=l.d0 

Jinv(2,l)=l.d0/b 

Jinv(2,2)=2.d0*yin/b 

return 

and 

Subroutine  »ap(xin,yin,xdin,ydin,xout ,yout , xdout , ydout ,a,b) 

Double  Precision  xin.yin.xout ,yout  ,a,b 

xout=a-xin*xin+b*yin 

yout-xin 

return 

End 
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